{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 Wine Classification with Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.special import *\n",
    "from save_csv import *\n",
    "\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shuffle_train_val_split(val_amt=0, percent=0):\n",
    "    data = loadmat(\"data.mat\")\n",
    "    total_num_ex = data[\"X\"].shape[0]\n",
    "    shuffle = np.random.permutation(total_num_ex) #shuffle before split\n",
    "    data_xtrain, data_ytrain = data[\"X\"][shuffle], \\\n",
    "    data[\"y\"][shuffle] #split\n",
    "    if not val_amt: #split by fixed number or percentage\n",
    "        val_amt = int(percent * total_num_ex)\n",
    "    data_xval, data_yval = data_xtrain[:val_amt,:], data_ytrain[:val_amt,:]\n",
    "    data_xtrain, data_ytrain = data_xtrain[val_amt:,:], data_ytrain[val_amt:,:]\n",
    "    return data_xtrain, data_ytrain, data_xval, data_yval\n",
    "\n",
    "def sigmoid(values):\n",
    "    return expit(values)\n",
    "\n",
    "def predict(X,w,mode=\"train\"):\n",
    "    predictions = sigmoid(np.dot(X,w))\n",
    "    if mode != \"train\": #validation or test\n",
    "        predictions = (predictions >= 0.5).astype(int)\n",
    "        \n",
    "    if mode == \"test\":\n",
    "        results_to_csv(predictions)\n",
    "        print(\"test predictions saved\")\n",
    "    else:\n",
    "        return predictions\n",
    "\n",
    "def logres_cost_reg(X,y,theta,lambd):\n",
    "    m = X.shape[0]\n",
    "    regularization = (lambd/(2*m)) * np.sum(np.square(theta[1:]))\n",
    "    predictions = predict(X,theta)\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    J = ((-1.0/m) * (np.sum(np.dot(-y, np.log(predictions)) + \\\n",
    "                xlog1py(1-y, -predictions)))) + regularization\n",
    "    return J\n",
    "\n",
    "def gradientdescent_reg_optimize(iterations,theta,lambd,alpha,GD_type=\"batch\",mode=\"train\"):\n",
    "    x, y = None, None\n",
    "    if mode == \"train\":\n",
    "        x, y = wine_xtrain, wine_ytrain\n",
    "    elif mode == \"test\":\n",
    "        x, y = wine_xfull, wine_yfull\n",
    "        \n",
    "    costs = [0] * iterations\n",
    "    if len(lambd) == 1:\n",
    "        lambd = [lambd[0]] * iterations\n",
    "    next_theta = theta.copy()\n",
    "    \n",
    "    for i in range(iterations):\n",
    "        ###\n",
    "        if GD_type == \"batch\":\n",
    "            new_theta = batchGD_reg_iteration(x,y,theta,lambd[i],alpha)\n",
    "        elif GD_type == \"stochastic\":\n",
    "            new_theta = stochasticGD_reg_iteration(x,y,theta,lambd[i],alpha)\n",
    "        costs[i] = logres_cost_reg(x,y,theta,lambd[i])\n",
    "        theta = new_theta\n",
    "        \n",
    "    return theta, costs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {},
   "outputs": [],
   "source": [
    "wine_xtrain, wine_ytrain, wine_xval, wine_yval = \\\n",
    "shuffle_train_val_split(percent=0.2)\n",
    "m1, n = wine_xtrain.shape\n",
    "m2, _ = wine_xval.shape\n",
    "wine_xtrain = np.concatenate([np.ones((m1, 1)), wine_xtrain], axis=1)\n",
    "wine_xval = np.concatenate([np.ones((m2, 1)), wine_xval], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 batch gradient descent update equation for logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![4.1 derivation](41updated.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.001\n",
    "lambd = [0.1] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "iterations = 1000\n",
    "\n",
    "def batchGD_reg_iteration(X,y,theta,lambd,alpha):\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    gradients = \\\n",
    "    (-1/len(y)*np.dot(X.T,(y-sigmoid(np.dot(X,theta))))).reshape(theta.shape)\n",
    "    gradients[1:] += (lambd/len(y) * theta[1:])\n",
    "    new_theta = theta - alpha*gradients\n",
    "    return new_theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEWCAYAAAC9qEq5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX9//HXGwJhC4R935RFUSliWLRaFUHRnxa0WrcWbFVa61Jb27r1W7W1X6W2dWn9alFxV7Su1I2K+4IgIIJsElkDCISwhSXr5/fHPdExJhAgyZ0kn+fjMY+Ze865M587meSTc+6Ze2RmOOecc9WhXtwBOOecqzs86TjnnKs2nnScc85VG086zjnnqo0nHeecc9XGk45zzrlq40nHJSVJyyUNr4bXuVHSY1X9OhUl6SFJN4fHx0haHGMsJqlXXK/vaidPOq7WkfS2pIvijmN/mdl7Zta3Mp6rupJ4dZB0gaT3447D7RtPOs5VEUkpccfgXLLxpOOS2SBJCyRtkvSgpEYAklpKeknShlD3kqQuoe7PwDHAPyXlSvpnKD9E0uuSciStk3Rdwus0lPSIpG2S5kvKKC8gSSdKWixpi6T/k/ROSa8q/Af+gaTbJeUAN0o6UNKbkjZKypb0uKT0hOc7XNLs8NpPAY0S6o6TlJWw3UnSs+G4l0m6IqHuRklPl3Uckh4FugH/Ce/J78o5tt9KWitpjaSflqpLlfRXSSvD+3evpMahrk34GWwO7+97kuqFuq6Sngsxbyz5eYS6n0paGH6GUyR1T6gzST+XtCTU363IwcC9wJHhWDaX97NyScrM/Oa3pLsBy4HPgK5AK+AD4OZQ1xr4AdAESAP+DbyQsO/bwEUJ22nAWuAqoj/qacCQUHcjsAs4BagP3AJ8VE5MbYCtwBlACvBLoKDktYALgELg8lDfGOgFjABSgbbAu8AdoX1DYAXwK6ABcGZ4vpLjPA7ICo/rAbOAP4T9DgCWAidV5DjC+zl8N+/3SGAdcCjQFHgCMKBXqL8DmBx+FmnAf4BbQt0tRImgQbgdAyjE8Slwe3jORsDRYZ/RQCZwcHivfg98mBCPAS8B6UQJcwMwMuF9fj/uz6jf9vF3O+4A/Oa3sm7hj+TPE7ZPAb4op+0AYFPCdumkcy7wSTn73ghMTdjuB+wsp+0YYFrCtoBVpZLOyj0c1+iSWIDvAWsAJdR/WE7SGVL6uYFrgQcrchwVSDoTgVsTtvuUJJ1wnNuBAxPqjwSWhcd/BF4sSVCl2mwAUsp4vVeBCxO26wE7gO5h20oSVNh+Grgm4X32pFNDbz7m7JLZqoTHK4BOAJKaEP33PBJoGerTJNU3s6Iynqcr8MVuXufLhMc7gEaSUsyssFS7TokxmZklDn+VETOS2gF3Ef33n0b0x3VTwvOttvCXNOE4y9Id6FRqOKk+8N4+HEdZOhH1pMqKoy1Rr3KWpJKykp4MwG1ESe+/oX6Cmd1K9L6vKOf1uwN3SvpbQpmAzgmvXfp4mlXgOFyS83M6Lpl1TXjcjahXANEwWV+iIbLmRD0GiP5oQfRfcqJVwIGVEM9aoEvJhqK/sF1KtSn92reEsv4h1h8lxLkW6KyEv+REx1mWVUQ9i/SEW5qZnVLB2Pd0Ofm1fPv9LpEN7AQOSXjtFmbWDMDMtpnZVWZ2AHAa8GtJJ4SYu5UzoWIV8LNSx9PYzD6shGNxScyTjktml0rqIqkVcB3wVChPI/ojuDnU3VBqv3VE5zxKvAR0kHRlOCGeJmnIPsTzMnCYpNHhD+mlQIc97JMG5IZYOwO/TaibRnQO6ApJKZLOAAaX8zwzgK2SrpbUWFJ9SYdKGlTB2Eu/J6U9DVwgqV/oSX71nppZMXAfcHvouSGps6STwuNTJfUKyXMrUBRuM4iS2a2SmkpqJOm74WnvBa6VdEh4jhaSztqLY+kiqWEF27sk4knHJbMngP8SnTBfCtwcyu8gOkmfDXwEvFZqvzuBM8Osp7vMbBvRyfzTiIZslgDH720wZpYNnAX8BdhIdN5kJpC3m91uAgYCW4iS1nMJz5dPNCnhAqIht7MT60u9dlGIfwCwjOjY7wdaVDD8W4Dfhxlmvynj+V8lel/fJDrB/2apJleH8o8kbQWmEvU2AXqH7VyiRPp/ZvZ2Qsy9gJVAVjhGzOx5YDwwKTzfZ8DJFTyWN4H5wJeSsiu4j0sS+uZwsnOuosK04CzgfDN7K+54nKsJvKfj3F6QdJKkdEmpREN+IuptOecqwJOOc3vnSKKZcNlEQ0ejzWxnvCE5V3P48Jpzzrlq4z0d55xz1ca/HFpKmzZtrEePHnGH4ZxzNcqsWbOyzaztntp50imlR48ezJw5M+4wnHOuRpFU3tU0vsGH15xzzlUbTzrOOeeqjScd55xz1caTjnPOuWrjScc551y18aTjnHOu2njScc45V238ezrOOVeLFRcb2/ML2Z5XRG5eITvyC8nNi7a35xWGukJy84oYfnA7+ndJr9J4POk451wSMzO27ixk8858tu4sZOuuArbuLGDrrgK27CwoVVYYygpCWSE7C8pawb1s7dJSPek451xtUlRsbNyeR/a2fHK255OzI5+c3DxydhSQsz2PTdsL2PjVfT6bd+RTWFz+hZnrCZo3bkDzRg1o3jiF5o0acGDbZjRvnEJaowY0S02haWp9mqamRI8bptAktX4oj8qaNKxPk4Yp1K+ncl+nsnjScc65SlBYVMyG3DzWb81j3dZdrN+Wx/pteWzYtisqC/fZuXmUlUMkSG/cgFZNG9KqaUN6tGnCwO7ptGwSbac3aUiLxg1o3iglSjKNG9CicQOaNqxPtFJ4zeBJxznnKiA3r5A1m3eyetNOVm/eGT0uud+0ky+37vpWMpGgddNU2qWl0q55Kv06Nqd980a0S0ulTbNUWjVtSOtmDWnZJEoq1dHTiJsnHeecIzrhvn5bHss3bmfFxu0s37gjus/ewerNO9mys+Ab7VPqiQ4tGtE5vTFDD2hN55aN6dCiEe3TGtGueSrt0hrRpllDUur7JOFESZd0JN1GtCJjPtEKjT8xs82h7lrgQqAIuMLMpoTykcCdQH3gfjO7NZT3BCYBrYDZwI/NLL96j8g5l0y27Cggc8M2lqzLZVn29pBkdrB843Z2FRR/1S6lnujWqgndWjfhiO4t6ZTemM4tG9M5vRGd05vQNi21TvRMKlvSJR3gdeBaMyuUNB64FrhaUj/gHOAQoBMwVVKfsM/dwAggC/hY0mQzWwCMB243s0mS7iVKWPdU8/E452KQsz2fJeu2sWR9Lpnrc1myPko067flfdWmYUo9urdqQvfWTTm6Vxu6t2lKj9ZN6NG6KR1bNPJeShVIuqRjZv9N2PwIODM8HgVMMrM8YJmkTGBwqMs0s6UAkiYBoyQtBIYB54U2DwM34knHuVqlqNhYvnE789dsZcGarSxYG91n536dXJqlptCrXTOO7dOWXu2a0bt9M3q3S6NzemPqeW+lWiVd0inlp8BT4XFnoiRUIiuUAawqVT4EaA1sNrPCMtp/g6RxwDiAbt26VUrgzrnKV1hUzOJ12/h01RYWrN3C/DVbWbR221ffRWlQX/Rpn8bxfdvSt0Mafdqn0bt9Mzo0b1SjZnjVZrEkHUlTgQ5lVF1vZi+GNtcDhcDjJbuV0d4o+1I+tpv23y40mwBMAMjIyCh/QrxzrtqYGWu37GLOqs3RbeVm5q3e8lWCSWuUQr+OzTlncFcO6dSCfh2b06tdMxqm+JBYMosl6ZjZ8N3VSxoLnAqcYGYlSSAL6JrQrAuwJjwuqzwbSJeUEno7ie2dc0mmsKiYhWu3MX3ZRj5ensMnKzd/df6lYf16HNI5SjADuqYzoGs63Vo18d5LDZR0w2thJtrVwLFmtiOhajLwhKS/E00k6A3MIOrR9A4z1VYTTTY4z8xM0ltE54QmAWOBF6vvSJxzu5NfWMy81ZuZviyH6UtzmLViE7l50Wh499ZNOOrA1lGC6daSgzumkZpSP+aIXWVIuqQD/BNIBV4P/8V8ZGY/N7P5kp4GFhANu11qZkUAki4DphBNmZ5oZvPDc10NTJJ0M/AJ8ED1HopzrkRxsbFg7VbeW5LN+5kbmLVi01dTlHu3a8bowzsxuGdrBvdoRYcWjWKO1lUVfT165SA6pzNz5sy4w3CuVli3dRfvLcnmvSUbeH9JNhu3R1+TO6hDGkce2JohPVsxqEcrWjdLjTlSt78kzTKzjD21S8aejnOuhioqNuas2sTrC9bz1qL1LF63DYA2zRryvT5tOaZ3G47u1YZ2zb0nU1d50nHO7Zcd+YW8tySbqQvW8eai9Wzcnk9KPTG4ZyuuGXgQx/Ruw8Edmvv3YRzgScc5tw+27ChgyvwvmTL/S97PzCavsJi0Rikc37cdw/u157i+bWneqEHcYbok5EnHOVch23YV8PqCdbw0dy3vLdlAQZHROb0x5w7uxoh+7RncsxUN/LIxbg886TjnyrUjv5CpC9fz0qdrePvzDeQXFtOpRSMuOKoHp/bvRP8uLfy7Mm6veNJxzn1DcbHx8fIcnpmVxSvz1rI9v4j2zVM5f0g3Tu3ficO7pvv5GbfPPOk45wBYlbODZ2dn8ezsLFbl7KRpw/qcclhHzhjYhSE9W3micZXCk45zdVheYRGvzvuSJ2esZPqyHCQ46sDW/Gp4H0Ye2oEmDf1PhKtc/olyrg5anr2dJ2es5N+zssjZnk/31k34zYl9OH1gFzqnN447PFeLedJxro4oLCpm6sL1PD59Be8tyaZ+PTHi4Pb8aGh3jjqwtQ+fuWrhSce5Wi5nez6Pf7SCx6avYN3WPDq2aMSvR/Th7EFdae9XBnDVzJOOc7XUknXbmPjBMp6bvZq8wmK+16ctfx7dneP6tvVlmF1sPOk4V4uYGe8uyeaB95fx7ucbSE2pxxkDu/DT7/agd/u0uMNzzpOOc7VBQVExL3yymgnvLmXJ+lzapqVy1Yg+nD+0O62aNow7POe+4knHuRpsV0ERT328ignvLmX15p0c3LE5fzvrO5z6nY6+6JlLSp50nKuBtu4q4LGPVjDx/WVk5+aT0b0lN59+KMf1aeuXpXFJLWmTjqTfALcBbc0sW9Fv0p3AKcAO4AIzmx3ajgV+H3a92cweDuVHAA8BjYFXgF+ar1rnarCc7flMfH8ZD09bzrZdhRzbpy2XHt+LwT1bxR2acxWSlElHUldgBLAyofhkoHe4DQHuAYZIagXcAGQABsySNNnMNoU244CPiJLOSODV6joO5yrLlh0FTHjvCx76YDk7Coo4+dAO/OK4XhzauUXcoTm3V5Iy6QC3A78DXkwoGwU8EnoqH0lKl9QROA543cxyACS9DoyU9DbQ3MymhfJHgNF40nE1yLZdBUx8fzn3v7+UbbsKObV/R64c3pte7XwmmquZki7pSPo+sNrMPi01Nt0ZWJWwnRXKdleeVUZ5Wa85jqhHRLdu3fbzCJzbfzvyC3now+VMeHcpm3cUcGK/9vxqRB8O7tg87tCc2y+xJB1JU4EOZVRdD1wHnFjWbmWU2T6Uf7vQbAIwASAjI8PP+bjY5BcW88T0FfzzrUyyc/M5vm9bfj2iL4d18WE0VzvEknTMbHhZ5ZIOA3oCJb2cLsBsSYOJeipdE5p3AdaE8uNKlb8dyruU0d65pGNmvPrZl/zltUUs37iDIw9ozb9+3JcjureMOzTnKlVSDa+Z2TygXcm2pOVARpi9Nhm4TNIkookEW8xsraQpwP9KKvntPBG41sxyJG2TNBSYDowB/lGdx+NcRcxakcOfX17I7JWb6dO+GQ9eMIjj+vrUZ1c7JVXS2YNXiKZLZxJNmf4JQEgufwI+Du3+WDKpALiEr6dMv4pPInBJZFn2dsa/uojX5n9Ju7RUxv/gMH4wsItfF83VavKvrXxTRkaGzZw5M+4wXC22dVcBd01dwkMfLic1pR4/O/ZALjqmpy+Y5mo0SbPMLGNP7fxT7lw1KS42np2dxfjXFrNxex5nZ3TlqhP70jYtNe7QnKs2nnScqwafrtrMDZPnM2fVZg7vls7ECzLo3yU97rCcq3aedJyrQhu25XHblEU8PTOLNs1S+dtZ3+H0wzv7Kp2uzvKk41wVKC42npixkvGvLWJnfhHjvncAlw/rRVqjBnGH5lysPOk4V8kWf7mN656fx6wVmzjygNb8afSh9GrXLO6wnEsKnnScqyS7Cor4x5tL+Nc7S2nWKIXbzuzPmUd08e/bOJfAk45zleCDzGyuf34eyzfu4IyBnbn+lINp3cxnpTlXmicd5/bDlh0F/PGlBTw7O4serZvw+EVD+G6vNnGH5VzS8qTj3D56Y+E6rn1uHjnb87ns+F5cNqwXjRr4EtHO7Y4nHef20pYdBdz00nyem72agzqkMfGCQb6YmnMV5EnHub3w5qKod5Odm88Vw3px2bDeNEzxa6U5V1GedJyrgK27CrhpcnTupm/7NO4fM8jXuHFuH3jScW4PPl6ew5WT5vDl1l1cdnwvLj+hF6kpfu7GuX3hSce5chQUFXPXG0u4+61MurRswr9/fiQDu/mias7tD086zpVhxcbt/HLSHOas2syZR3Thxu8fQrNU/3Vxbn8l5RlQSZdLWixpvqS/JJRfKykz1J2UUD4ylGVKuiahvKek6ZKWSHpKUsPqPhZXs5gZz8zK4pQ732Pphlz+ed7h/PWs73jCca6SJN1vkqTjgVFAfzPLk9QulPcDzgEOAToBUyX1CbvdDYwAsoCPJU02swXAeOB2M5sk6V7gQuCe6j0iV1Nszyvkuufn8eKcNQzp2Yrbzx5Ap/TGcYflXK2SdEmHaInpW80sD8DM1ofyUcCkUL5MUiYwONRlmtlSAEmTgFGSFgLDgPNCm4eBG/Gk48qw+Mtt/OLxWSzL3s5VI/rwi+N7Ud+XH3Cu0iXj8Fof4JgwLPaOpEGhvDOwKqFdVigrr7w1sNnMCkuVf4ukcZJmSpq5YcOGSjwUVxM8MyuLUXe/z5adhTx20RAuP6G3JxznqkgsPR1JU4EOZVRdTxRTS2AoMAh4WtIBQFl/BYyyE6ftpv23C80mABMAMjIyymzjap+d+UXcMPkznp6ZxdADWnHXuYfTLq1R3GE5V6vFknTMbHh5dZIuAZ4zMwNmSCoG2hD1VLomNO0CrAmPyyrPBtIlpYTeTmJ7V8ct3ZDLLx6fzaIvt3H5sF5cObyP926cqwbJOLz2AtG5GMJEgYZECWQycI6kVEk9gd7ADOBjoHeYqdaQaLLB5JC03gLODM87FnixWo/EJaU3Fq5j1D8/YN3WXTz0k0FcdWJfTzjOVZNknEgwEZgo6TMgHxgbEsh8SU8DC4BC4FIzKwKQdBkwBagPTDSz+eG5rgYmSboZ+AR4oHoPxSWT4mLj7rcy+fvUz+nXsTn/+vERdGnZJO6wnKtTFP09dyUyMjJs5syZcYfhKlluXiG/efpTXpv/JaMHdOKWM/rTuKFfysa5yiJplpll7KldMvZ0nKtUy7O3M+7RmWSuz+X3/+9gLjy6py8h7VxMPOm4Wu2dzzdw+ROzqVdPPPLTIRzd21f1dC5OnnRcrfXItOXcOHk+fdqncd+YDLq28vM3zsXNk46rdQqLirn55YU89OFyhh/cjjvPOZymfu0055KC/ya6WiU3r5DLn5jNW4s3cOHRPbnulIN9OrRzScSTjqs1Vm/eyYUPfcyS9bncPPpQfjS0e9whOedK8aTjaoU5qzZz0cMzySso4qGfDOKY3m3jDsk5VwZPOq7Ge2PhOi59YjZtmqXy5MVD6N0+Le6QnHPl8KTjarRJM1Zy3fPzOLRzCx4YO4i2aalxh+Sc2w1POq5GMjP+8WYmf3/9c77Xpy33nD/QZ6g5VwP4b6mrcYqKjf958TOemL6SMwZ2ZvwP+tOgfjJeu9Y5V5onHVej7Coo4oonP+G/C9ZxyXEH8ruT+volbZyrQTzpuBpjy84CLnzoY2at3MSNp/Xjgu/2jDsk59xe8qTjaoTs3Dx+/MAMMtdv4x/nHs6p/TvFHZJzbh940nFJb+2WnZx//3TWbN7J/WMHcWwf/w6OczWVJx2X1FZs3M55901n684CHvnpEAb3bBV3SM65/ZB0U34kDZD0kaQ5kmZKGhzKJekuSZmS5koamLDPWElLwm1sQvkRkuaFfe6Sn3GuUT5ft42z7p3GjvxCnrh4qCcc52qBpEs6wF+Am8xsAPCHsA1wMtA73MYB9wBIagXcAAwBBgM3SGoZ9rkntC3Zb2Q1HYPbT3OzNvPDf00D4KmfHclhXVrEHJFzrjJUeHhN0lFAj8R9zOyRKojJgObhcQtgTXg8CnjEovW1P5KULqkjcBzwupnlhDhfB0ZKehtobmbTQvkjwGjg1SqI2VWi2Ss3MfaBGbRo0oDHLxpC99ZN4w7JOVdJKpR0JD0KHAjMAYpCsQFVkXSuBKZI+itRT+yoUN4ZWJXQLiuU7a48q4zyb5E0jqhHRLdu3fb/CNw+m71yE2MemEHrZg158uKhdEpvHHdIzrlKVNGeTgbQL/Qy9pukqUCHMqquB04AfmVmz0r6IfAAMBwo63yM7UP5twvNJgATADIyMirlGN3em7ViE2MnRgln0rihdGzhCce52qaiSeczoiSxtjJe1MyGl1cXhsF+GTb/DdwfHmcBXROadiEaessiGmJLLH87lHcpo71LQiUJp02zhjzpCce5WquiEwnaAAskTZE0ueRWRTGtAY4Nj4cBS8LjycCYMIttKLDFzNYCU4ATJbUMEwhOBKaEum2ShoZZa2OAF6soZrcfZq3I+SrhTBp3pCcc52qxivZ0bqzKIEq5GLhTUgqwi3CuBXgFOAXIBHYAPwEwsxxJfwI+Du3+WDKpALgEeAhoTDSBwCcRJJnZKzcxduLHtE1L5cmLh9KhRaO4Q3LOVSFV9DSNpPbAoLA5w8zWV1lUMcrIyLCZM2fGHUadMH/NFs6Z8BGtm0Y9HE84ztVckmaZWcae2lVoeC2c0J8BnAX8EJgu6cz9C9HVZZnrcxnzwAzSUlN47KIhnnCcqyMqOrx2PTCopHcjqS0wFXimqgJztdeqnB386P7pSOKxi4bQpWWTuENyzlWTik4kqFdqOG3jXuzr3FfWbd3F+fdPZ2dBEY9eOJgD2jaLOyTnXDWqaE/nNUlTgCfD9tlEJ/adq7Cc7fn86P7pbMzN47GLhnBwx+Z73sk5V6tUKOmY2W8l/QD4LtGXLieY2fNVGpmrVXLzChk7cQYrc3bw0E8Gc3i3lnveyTlX61T42mtm9izwbBXG4mqp/MJiLnlsFgvWbuW+MUdw5IGt4w7JOReT3SYdSe+b2dGStvHNS8gIMDPz8RG3W2bGNc/O5b0l2fzlzP4MO6h93CE552K026RjZkeH+7TqCcfVNuNfW8xzn6zmqhF9+GFG1z3v4Jyr1Sr6PZ1HK1LmXKKHPljGve98wflDunHZsF5xh+OcSwIVnfZ8SOJGuETNEZUfjqstXp67lpteWsCJ/drzx1GH4ou2OudgD0lH0rXhfE5/SVvDbRuwDr94pivH9KUb+dVTcziiW0vuOvdw6tfzhOOci+w26ZjZLeF8zm1m1jzc0systZldW00xuhpkWfZ2fvbYLLq0asz9YzNo1KB+3CE555JIRb+nc21YNqA30Cih/N2qCszVPJt35HPhQx8j4MELBpHepGHcITnnkkxFl6u+iGhhtS5ES1YPBaYRrXfjHPmFxfz8sVlkbdrJ4xcPoXvrpnGH5JxLQhWdSPBLomUNVpjZ8cDhwIYqi8rVKGbG71+Yx0dLcxh/5mEM6tEq7pCcc0mqoklnl5ntApCUamaLgL77+qKSzpI0X1KxpIxSdddKypS0WNJJCeUjQ1mmpGsSyntKmi5piaSnJDUsiTNsZ4b6Hvsar9u9e99ZytMzs7hiWC9OP7zLnndwztVZFU06WZLSgReA1yW9SLSs9L76DDgD+MY5IUn9gHOIpmiPBP5PUn1J9YG7gZOBfsC5oS3AeOB2M+sNbAIuDOUXApvMrBdwe2jnKtlrn61l/GuLOO07nfjViD5xh+OcS3IVSjpmdrqZbTazG4H/AR4ARu/ri5rZQjNbXEbVKGCSmeWZ2TKipakHh1ummS01s3xgEjBK0Zc/hvH1uj4PJ8Q1KmwT6k+Qf1mkUs3N2syVT81hYLd0bjuzv38Xxzm3R3tMOpLqSfqsZNvM3jGzyeGPf2XrDKxK2M4KZeWVtwY2m1lhqfJvPFeo3xLau0qwftsuxj0yi9ZNU5kwxqdGO+cqZo+z18ysWNKnkrqZ2cqKPrGkqUCHMqquN7Pyvlha1r/KRtnJ0XbTfnfP9e0XlcYB4wC6detWTmiuRF5hEZc8NpstOwt49pKjaNMsNe6QnHM1REWXNugIzJc0A9heUmhm3y9vBzMbvg/xZAGJV4XswtfnjsoqzwbSJaWE3kxi+5LnygqX7WkB5JQT6wRgAkBGRkaZiclFzIwbXpzPrBWbuPu8gfTr5Bcad85VXEWTzk1VGsXXJgNPSPo70Inoy6gziHotvSX1BFYTTTY4z8xM0lvAmUTnecby9eV5JoftaaH+TTPzhLKfHvtoBZM+XsWlxx/I/+vfMe5wnHM1TEWvSPCOpO5AbzObKqkJsM+D+JJOB/4BtAVeljTHzE4ys/mSngYWAIXApWZWFPa5DJgSXneimc0PT3c1MEnSzcAnRJMcCPePSsok6uGcs6/xusj0pRu56T8LGHZQO64asc8z5p1zdZgq8s+/pIuJznm0MrMDJfUG7jWzE6o6wOqWkZFhM2fOjDuMpLN6806+/4/3adGkAS9c+l2aN2oQd0jOuSQiaZaZZeypXUW/p3Mp8F1gK4CZLQHa7Xt4ribZmV/EuEdmkl9YzH1jMjzhOOf2WUWTTl7iFOlwYt7Pj9QBZsb1L8xjwdqt3HnuAA5s2yzukJxzNVhFk847kq4DGksaAfwb+E/VheWSxZMzVvHc7NX88oTeDDuofdzhOOdquIomnWuILvA5D/gZ8IqZXV9lUbmkMDdrMzdOns/3+rTlimG94w7HOVcLVHTK9OVmdidwX0mBpF+GMlcLbd6RzyWPzaZtWip3nD2Aer76p3MpLo07AAATyUlEQVSuElS0pzO2jLILKjEOl0SKi41fPTWH9dt2cff5A2nV1Bdjc85Vjt32dCSdC5wH9JQ0OaEqDdhYlYG5+Nz9ViZvLd7An0YdwoCu6XGH45yrRfY0vPYhsBZoA/wtoXwbMLeqgnLxeX9JNn+f+jmjBnTiR0O7xx2Oc66W2W3SMbMVwArgyOoJx8Xpyy27uGLSJ/Rq24xbzjjMlypwzlW6PQ2vvW9mR0vaxje/lyPAzMyv9lhLFBUbV0z6hF0FRdzzoyNo0rCic0ycc67i9tTTOTrcp1VPOC4u/3hzCTOW5fC3s75Dr3b+BVDnXNWo6Ow1V4t9tHQjd72xhDMO78wPjugSdzjOuVrMk04dl7M9nysnzaF766b8cfShcYfjnKvlfOC+DjMzfvvvT8nZns9zY4+iWap/HJxzVct7OnXYxA+W88ai9Vx3ykEc2rlF3OE45+oATzp11LysLdz66kKGH9yesUf1iDsc51wdEUvSkXSWpPmSiiVlJJSPkDRL0rxwPyyh7ohQninpLoUvkUhqJel1SUvCfctQrtAuU9JcSQOr/0iTU25eIZc/OZs2zVK57cz+/n0c51y1iaun8xlwBvBuqfJs4DQzO4zoem+PJtTdQ7R6ae9wGxnKrwHeMLPewBthG+DkhLbjwv4O+NN/FrAiZwd3nD2Aln5dNedcNYol6ZjZQjNbXEb5J2a2JmzOBxpJSpXUEWhuZtMsWl/7EWB0aDcKeDg8frhU+SMW+QhID89Tp02Z/yVPzVzFJcceyJADWscdjnOujknmczo/AD4xszygM5CVUJcVygDam9lagHBfsox2Z2BVOft8g6RxkmZKmrlhw4ZKPITksn7rLq55di6Hdm7OlcP7xB2Oc64OqrI5spKmAh3KqLrezF7cw76HAOOBE0uKymi2p+WyK7yPmU0AJgBkZGTUymW4zYzfPTuXHflF3HH2ABqmJPP/G8652qrKko6ZDd+X/SR1AZ4HxpjZF6E4C0j8qnwXoGQYbp2kjma2NgyfrU/Yp2s5+9Q5j320grcXb+CPow6hVzu/qpFzLh5J9e+upHTgZeBaM/ugpDwMm22TNDTMWhsDlPSWJvP1InNjS5WPCbPYhgJbSobh6prM9du4+eWFHNunLT/25QqcczGKa8r06ZKyiJZMeFnSlFB1GdAL+B9Jc8Kt5BzNJcD9QCbwBfBqKL8VGCFpCTAibAO8AiwN7e8DflHFh5WU8guLufKpOTRpWN+nRzvnYqdoMpgrkZGRYTNnzow7jEpz25RF3P3WF9z7oyMYeWhZp9icc27/SZplZhl7apdUw2uucs1akcM9b3/BDzO6eMJxziUFTzq11M78In7z77l0bNGYP5x2SNzhOOcc4FeZrrVum7KYZdnbeeKiIX71aOdc0vCeTi00felGHvxwGWOO7M5RvdrEHY5zzn3Fk04tsyO/kN8+M5euLZtw9ciD4g7HOee+wcddapnxry5i1aYdTLp4KE19WM05l2S8p1OLfPhFNg9PW8FPjurpF/N0ziUlTzq1RG5eIb97Zi492zTltyf1jTsc55wrk4+/1BK3vLKQ1Zt38szPj6Rxw/pxh+Occ2Xynk4t8N6SDTw+fSUXH3MAR3RvFXc4zjlXLk86NVxuXiHXPDuPA9o25dcjfI0c51xy8+G1Gu6vUxazZstO/v2zI2nUwIfVnHPJzXs6NdisFZt4eNpyxgztTkYPH1ZzziU/Tzo1VF5hEVc/O5eOzRvxW/8SqHOuhvDhtRrq7re+IHN9Lg/+ZJBfW805V2PEtYjbWZLmSyqW9K31FyR1k5Qr6TcJZSMlLZaUKemahPKekqZLWiLpKUkNQ3lq2M4M9T2q49iqw6Ivt/J/b2Vy+uGdOb5vuz3v4JxzSSKu4bXPgDOAd8upv52vVwZFUn3gbuBkoB9wrqR+oXo8cLuZ9QY2AReG8guBTWbWKzzf+Mo+iDgUFRtXPzuP5o0b8D+n9tvzDs45l0RiSTpmttDMFpdVJ2k00TLT8xOKBwOZZrbUzPKBScAoRWsvDwOeCe0eBkaHx6PCNqH+BNWCtZof/GAZn67azI3fP4RWTRvGHY5zzu2VpJpIIKkpcDVwU6mqzsCqhO2sUNYa2GxmhaXKv7FPqN8S2tdYKzfu4K//XcwJB7XjtP4d4w7HOef2WpWdgZY0FShrjeTrzezFcna7iWioLLdUp6SsHortpnx3+5QV6zhgHEC3bt3KCS1eZsZ1z88jpV49bj79UGpBp805VwdVWdIxs+H7sNsQ4ExJfwHSgWJJu4BZQNeEdl2ANUA2kC4pJfRmSsoh6vV0BbIkpQAtgJxyYp0ATADIyMgoMzHF7ZlZWbyfmc2fRh9KxxaN4w7HOef2SVLNtTWzY0oeS7oRyDWzf4ak0VtST2A1cA5wnpmZpLeAM4nO84wFSnpRk8P2tFD/ppklZULZk425efz5lYVkdG/J+YOTsyfmnHMVEdeU6dMlZQFHAi9LmrK79qEXcxkwBVgIPG1mJRMNrgZ+LSmT6JzNA6H8AaB1KP81cA011P++sojcXYX87xmHUa+eD6s552quWHo6ZvY88Pwe2txYavsV4JUy2i0lmt1WunwXcNZ+BZoEpn2xkWdnZ/GL4w6kT/u0uMNxzrn9klSz19w35RUWcf0L8+jWqgmXD+sddzjOObffkuqcjvume99eytIN23n4p4N9YTbnXK3gPZ0ktSx7O3e/ncmp/TtybJ+2cYfjnHOVwpNOEjIzfv/CPFLr1+MPfqkb51wt4kknCb0wZzUfZG7kdycfRLvmjeIOxznnKo0nnSSzeUc+N7+0kAFd0/07Oc65WscnEiSZ8a8tYvPOAh493b+T45yrfbynk0RmLs/hyRmruPDonvTr1DzucJxzrtJ50kkShUXF/P6Fz+ic3pgrh/t3cpxztZMnnSTxyLQVLPpyG/9zaj+aNPRRT+dc7eRJJwms37qL21//nGP7tOWkQ9rHHY5zzlUZTzpJ4JZXF5FXWMyN3z/E18lxztVqnnRiNn3pRp7/ZDU/O/YAerZpGnc4zjlXpTzpxKigqJg/vDifzumN+cVxveIOxznnqpwnnRg9Mm0Fi9dt4w+n9fMLejrn6gRPOjFZv3UXd7z+Ocf1bcuJ/XzygHOubohr5dCzJM2XVCwpo1Rdf0nTQv08SY1C+RFhO1PSXQpn3CW1kvS6pCXhvmUoV2iXKWmupIHVf6Tl+2rywGk+ecA5V3fE1dP5DDgDeDexUFIK8BjwczM7BDgOKAjV9wDjgN7hNjKUXwO8YWa9gTf4elnqkxPajgv7J4XEyQM9fPKAc64OiSXpmNlCM1tcRtWJwFwz+zS022hmRZI6As3NbJqZGfAIMDrsMwp4ODx+uFT5Ixb5CEgPzxMrnzzgnKvLku2cTh/AJE2RNFvS70J5ZyAroV1WKANob2ZrAcJ9u4R9VpWzzzdIGidppqSZGzZsqKRDKVvJ5IEbfPKAc64OqrLrrUiaCnQoo+p6M3txN/EcDQwCdgBvSJoFbC2jre0phIruY2YTgAkAGRkZe3refVZy5YHj+7ZlhE8ecM7VQVWWdMxs+D7slgW8Y2bZAJJeAQYSnefpktCuC7AmPF4nqaOZrQ3DZ+sTnqtrOfvE4tbXFpFfWMwNPnnAOVdHJdvw2hSgv6QmYVLBscCCMGy2TdLQMGttDFDSW5oMjA2Px5YqHxNmsQ0FtpQMw8Xhk5WbeG72ai46pqdPHnDO1VlxTZk+XVIWcCTwsqQpAGa2Cfg78DEwB5htZi+H3S4B7gcygS+AV0P5rcAISUuAEWEb4BVgaWh/H/CLqj6u8hQXGzf+ZwHt0lL5xfE+ecA5V3fFcg19M3seeL6cuseIhtNKl88EDi2jfCNwQhnlBly638FWguc+Wc2nqzbz9x9+h2apvmyBc67uSrbhtVonN6+Q8a8tYkDXdEYPKHPynHPO1Rn+b3cV++ebmWzYlsd9YzKoV88nDzjn6jbv6VSh5dnbmfj+Mn4wsAsDuqbHHY5zzsXOk04VuvnlhTSoL64e2TfuUJxzLil40qki736+gakL13HZsN60a94o7nCccy4peNKpAgVFxfzxpQV0b92Enx7dI+5wnHMuaXjSqQKPTltB5vpcfv//+pGa4tdXc865Ep50KtnG3Dxun/o5x/Ruw/CD2+15B+ecq0M86VSyv73+OTvyi7jhtH5+fTXnnCvFk04lmr9mC0/OWMmYI7vTq11a3OE451zS8aRTScyMm/6zgJZNGnLlCX3iDsc555KSJ51K8vK8tcxYlsNVJ/ahRZMGcYfjnHNJyZNOJWmWmsKIfu05Z1C3uENxzrmk5ddeqyTH9W3HcX19tppzzu2O93Scc85Vm7gWcTtL0nxJxZIyEsobSHpY0jxJCyVdm1A3UtJiSZmSrkko7ylpuqQlkp6S1DCUp4btzFDfozqP0Tnn3LfF1dP5DDgDeLdU+VlAqpkdBhwB/ExSD0n1gbuBk4F+wLmS+oV9xgO3m1lvYBNwYSi/ENhkZr2A20M755xzMYol6ZjZQjNbXFYV0FRSCtAYyAe2AoOBTDNbamb5wCRglKJvXw4Dngn7PwyMDo9HhW1C/Qnyb2s651ysku2czjPAdmAtsBL4q5nlAJ2BVQntskJZa2CzmRWWKidxn1C/JbR3zjkXkyqbvSZpKtChjKrrzezFcnYbDBQBnYCWwHvhecrqodhuytlDXelYxwHjALp18ynPzjlXVaos6ZjZ8H3Y7TzgNTMrANZL+gDIIOqxdE1o1wVYA2QD6ZJSQm+mpByiXk9XICsM17UAcsqJdQIwASAjI6PMxOScc27/Jdvw2kpgmCJNgaHAIuBjoHeYqdYQOAeYbGYGvAWcGfYfC5T0oiaHbUL9m6G9c865mCiOv8OSTgf+AbQFNgNzzOwkSc2AB4lmqAl40MxuC/ucAtwB1AcmmtmfQ/kBRBMLWgGfAD8yszxJjYBHgcOJejjnmNnSCsS2AVhRmcdbA7Qh6jXWZXX9Pajrxw/+HsD+vQfdzaztnhrFknRccpE008wy9tyy9qrr70FdP37w9wCq5z1ItuE155xztZgnHeecc9XGk46DMHOvjqvr70FdP37w9wCq4T3wczrOOeeqjfd0nHPOVRtPOs4556qNJ51aTlJXSW+FpSLmS/plKG8l6fWwJMTrklqGckm6KywJMVfSwHiPoPJIqi/pE0kvhe06tSyGpHRJz0haFD4PR9alz4GkX4Xfgc8kPSmpUV34DEiaKGm9pM8Syvb65y5pbGi/RNLYsl6rIjzp1H6FwFVmdjDRFR4uDctCXAO8EZaEeCNsQ7R8RO9wGwfcU/0hV5lfAgsTtuvashh3El1m6iDgO0TvRZ34HEjqDFwBZJjZoURfMj+HuvEZeAgYWapsr37ukloBNwBDiK6ReUNJotprZua3OnQjukzQCGAx0DGUdQQWh8f/As5NaP9Vu5p8I7ou3xtES2G8RHTFi2wgJdQfCUwJj6cAR4bHKaGd4j6G/Tz+5sCy0sdRVz4HfH3V+VbhZ/oScFJd+QwAPYDP9vXnDpwL/Cuh/Bvt9ubmPZ06JAwRHA5MB9qb2VqAcN8uNCtvGYma7g7gd0Bx2K5ry2IcAGwAHgxDjPeH6xvWic+Bma0G/kp0fce1RD/TWdStz0Civf25V9rnwZNOHRGua/cscKWZbd1d0zLKavS8ekmnAuvNbFZicRlN93pZjBokBRgI3GNmhxOtW3XNbtrXqvcgDAWNAnoSLZ3SlGgoqbTa/BmoiH1ZRmaveNKpAyQ1IEo4j5vZc6F4naSOob4jsD6UlywJUSJxuYia6rvA9yUtJ7o47DCink96WPYCyl4Wgz0ti1GDZAFZZjY9bD9DlITqyudgOLDMzDZYtHTKc8BR1K3PQKK9/blX2ufBk04tJ0nAA8BCM/t7QlXi0g+ll4QYE2axDAW2lHTDayozu9bMuphZD6KTx2+a2fnUoWUxzOxLYJWkvqHoBGABdedzsBIYKqlJ+J0oOf468xkoZW9/7lOAEyW1DL3GE0PZ3ov7BJffqvYGHE3UDZ4LzAm3U4jGp98AloT7VqG9gLuBL4B5RLN9Yj+OSnw/jgNeCo8PAGYAmcC/gdRQ3ihsZ4b6A+KOu5KOfQAwM3wWXiBanbfOfA6Am4jW5/qMaNmT1LrwGQCeJDqPVUDUY7lwX37uwE/D+5EJ/GRf4/HL4DjnnKs2PrzmnHOu2njScc45V2086TjnnKs2nnScc85VG086zjnnqo0nHeeSgKQOkiZJ+kLSAkmvSOqzl89xXVXF51xl8SnTzsUsfFnxQ+BhM7s3lA0A0szsvb14nlwza1ZFYTpXKbyn41z8jgcKShIOgJnNAd6XdFtY/2WepLMhumyJpHclzQl1x0i6FWgcyh6P6Tic26OUPTdxzlWxQ4mueFzaGURXEfgO0Ab4WNK7wHlEl+D/s6T6QBMze0/SZWY2oNqidm4feNJxLnkdDTxpZkVEF2h8BxgEfAxMDBdyfSH0ipyrEXx4zbn4zQeOKKO8rMvJY2bvAt8DVgOPShpThbE5V6k86TgXvzeBVEkXlxRIGkS0fPLZkupLakuUaGZI6k60PtB9RFcQL1nHviD0fpxLWj685lzMzMwknQ7cIekaYBewHLgSaAZ8SnSl8N+Z2ZeSxgK/lVQA5AIlPZ0JwFxJsy1ausG5pONTpp1zzlUbH15zzjlXbTzpOOecqzaedJxzzlUbTzrOOeeqjScd55xz1caTjnPOuWrjScc551y1+f8u4A2A2VPemAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta, costs = gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"batch\")\n",
    "plt.title(\"batch gradient descent\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 stochastic gradient descent update equation for logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 420,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "lambd = [0.01] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "iterations = 1000\n",
    "\n",
    "def stochasticGD_reg_iteration(X,y,theta,lambd,alpha):\n",
    "    new_theta = theta.copy()\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    np.random.shuffle(X)\n",
    "    for i in range(len(y)):\n",
    "        gradients = np.dot(y[i]-sigmoid(np.dot(X[i,:],theta)),X[i,:])\n",
    "        gradients[1:] += (lambd/len(y) * theta[1:])\n",
    "        new_theta = theta - alpha*gradients\n",
    "    return new_theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/cs189/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xb26931c88>]"
      ]
     },
     "execution_count": 429,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGV1JREFUeJzt3Xu0XGWd5vHvQ6IBBLkGgYQQNLTTwbZRDzhMYy+UW3CNQis2aM8yjBd61mi3l9E2gNNc1AVoO9i91HEiOk3jBVBbje2F5iKo7QVOEJWodMJNAgiBIIIICP7mj9rB4kydnMo5u07lmO9nrVq197vf2vv31knqOXu/dapSVUiSNFVbDbsASdLvBwNFktQKA0WS1AoDRZLUCgNFktQKA0WS1AoDRZu1JJVk0TQc56tJlg76OJsqyT8meXez/IIk1w+xlmn5WWjmMlDUuiSnJfnEsOsYT6/6quqoqjpvWDX1o6q+WVXPbGNfSW5Oclgb+xq2JCck+daw65CBIk2bJLOHXYM0SAaKJi3JO5LcluT+JNcnOTTJEuBk4LgkDyT5QdN3zyQrkqxPsibJ67v2MyvJyUluaPa1MsleXYc6LMnqJPcm+VCSNI97RpLLk9yT5O4kn0yy4yTruyLJ67oe+/okP2ke++Mkzx3nOTii2fd9ST6c5MoN+2l+c/63JOckWQ+c1kfNz0lyTXPcC4Gtu7YdkmRt1/qeST6XZF2Sm5L8dde205JclOSfmn2tSjLSbDsfWAB8qXkO/macsb09yR1Jbk/ymjHb5iT5uyQ/S3Jnko8k2abZtmuSf0nyi+bn/c0kWzXb9kryz03N9yT5YNc+X9M85/cmuTjJ3l3bKsl/G/vvIMkfAh8BDmrG8oteY9E0qSpv3jb5BjwTuBXYs1lfCDyjWT4N+MSY/lcCH6bzArk/sA44tNn2duBHzT4D/DGwS7OtgH8BdqTzIrgOWNJsWwQcDswB5gLfAD4wyfquAF7XLL8CuA04oKlnEbB3j+dgV+CXwMuA2cCbgN907ecE4FHgr5rt20xQ85OBW4C3AE8Cjm329+5m+yHA2mZ5K2Al8LfN454O3Agc2TXGh4AXA7OAM4HvdtV+M3DYRn6+S4A7gWcBTwE+1fwsFjXbPwCsAHYGtge+BJzZbDuTzov8k5rbC5rncRbwA+CcZp9bAwc3jzkGWAP8YfNcvRP4dlc9G/t3cALwrWH/n/BWBoq3yd2aF8a7gMOAJ43Z9oQXbGAv4DFg+662M4F/bJavB44e5zi14UWnWb8IWDZO32OA729qfU3bFfwuCC4G3tTHc/Bq4Dtd66ETYt2B8rMJ9tFd858CtwPp2v5tegfK88fuGzgJ+L9dY7y0a9ti4Ndd6zez8UD5OHBW1/ofND+LRc04f0UT0M32g4CbmuUzgC/ShM+YPuuA2T2O91XgtV3rWwEP0gT5xv4dGCibz81LXpqUqloDvJnOC9ddSS5Isuc43fcE1lfV/V1ttwDzmuW9gBs2crifdy0/CGwHkGS35ri3Jfkl8Ak6Zw2bWt9YE9WzwZ50AoTmmAWsHdPn1u6VjdXc7O+2Zj8b3DLOsfcG9mwuK/2iudRzMvC0rj5jn7et0/88zhPGNqaOucC2wMquY3+taQd4H52zjX9NcmOSZU37XsAtVfXoOOP5+679racTXPO6+vT8d6DNh4GiSauqT1XVwXReDAo4e8OmMV1vB3ZOsn1X2wI6l5Wg88L1jEmUcGZzrGdX1VOB/0LnRWhT6xur33ruAOZvWGnmduaP6TP2WBur+Q5g3oY5osaCjdR4U1Xt2HXbvqpe3Efdveoa6w46AdCrjruBXwP7dR17h6raDqCq7q+q/1FVTwdeArw1yaFNzQvGCbVbgb8cM55tqurbLYxF08RA0aQkeWaSFyWZQ+da/a/pXNaCzrX3hRsmYqvqVjqXbs5MsnWSZwOvBT7Z9D8XeFeSfZuJ1mcn2aWPMrYHHgB+kWQenbmYTa6vh3OBtyV5XlPPou4J4i5fBv4oyTHNi+QbgN0nWzPwHTpzLn+dZHaSlwEHjrOfq4BfpvPGg23SeWPDs5IcMMHxN7iTzrzLeC4CTkiyOMm2wKkbNlTVb4GPAuck2Q0gybwkRzbL/7l5zkJnjumx5nYVnaA6K8lTmn8Lf9Ls9iPASUn2a/axQ5JXbMJY5id5cp/9NSAGiiZrDnAWnd9Wfw7sRueSC8Bnmvt7klzTLL+SzsT47cDngVOr6pJm2/+i8wL2r3RegD5GZwJ7IqcDzwXuo/Pi/s9TqO9xVfUZ4D10JqLvB75AZ/J5bL+76Uzgvxe4h848xSjw8GRqrqpH6EzwnwDcCxw3Zkzdx36Mzm//+wM3NeM8F9hhI8fudibwzuYS09t67P+rdCbeL6dz+eryMV3e0bR/t7l0dymdN0IA7NusP0AnJD9cVVd01bwI+Bmdy4PHNcf7PJ0zyAua/V0HHNXnWC4HVgE/T3J3n4/RAOSJl2slTVZzxrMW+Iuq+vqw65Gmm2co0hQkOTLJjs2ltZPpzId8d8hlSUNhoEhTcxCdd4TdTedyzjFV9evhliQNh5e8JEmt8AxFktSKLerD6nbddddauHDhsMuQpBll5cqVd1fV3In6bVGBsnDhQkZHR4ddhiTNKEnG+8SGJ/CSlySpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFQaKJKkVBookqRUGiiSpFUMNlCRLklyfZE2SZT22z0lyYbP9e0kWjtm+IMkDSd42XTVLknobWqAkmQV8CDgKWAy8MsniMd1eC9xbVYuAc4Czx2w/B/jqoGuVJE1smGcoBwJrqurGqnoEuAA4ekyfo4HzmuXPAocmCUCSY4AbgVXTVK8kaSOGGSjzgFu71tc2bT37VNWjwH3ALkmeArwDOH2igyQ5McloktF169a1Urgk6f83zEBJj7bqs8/pwDlV9cBEB6mq5VU1UlUjc+fOnUSZkqR+zB7isdcCe3WtzwduH6fP2iSzgR2A9cDzgWOTvBfYEfhtkoeq6oODL1uS1MswA+VqYN8k+wC3AccDrxrTZwWwFPgOcCxweVUV8IINHZKcBjxgmEjScA0tUKrq0SRvBC4GZgEfr6pVSc4ARqtqBfAx4Pwka+icmRw/rHolSRuXzi/8W4aRkZEaHR0ddhmSNKMkWVlVIxP18y/lJUmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrRhqoCRZkuT6JGuSLOuxfU6SC5vt30uysGk/PMnKJD9q7l803bVLkp5oaIGSZBbwIeAoYDHwyiSLx3R7LXBvVS0CzgHObtrvBl5SVX8ELAXOn56qJUnjGeYZyoHAmqq6saoeAS4Ajh7T52jgvGb5s8ChSVJV36+q25v2VcDWSeZMS9WSpJ6GGSjzgFu71tc2bT37VNWjwH3ALmP6vBz4flU9PKA6JUl9mD3EY6dHW21KnyT70bkMdsS4B0lOBE4EWLBgwaZXKUnqyzDPUNYCe3WtzwduH69PktnADsD6Zn0+8Hng1VV1w3gHqarlVTVSVSNz585tsXxJUrdhBsrVwL5J9knyZOB4YMWYPivoTLoDHAtcXlWVZEfgy8BJVfVv01axJGlcQwuUZk7kjcDFwE+Ai6pqVZIzkry06fYxYJcka4C3AhveWvxGYBHwP5Nc29x2m+YhSJK6pGrstMXvr5GRkRodHR12GZI0oyRZWVUjE/XzL+UlSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa2Y3W/HJP8JWNj9mKr6pwHUJEmagfoKlCTnA88ArgUea5oLMFAkSUD/ZygjwOLakr4vWJK0SfqdQ7kO2H2QhUiSZrZ+z1B2BX6c5Crg4Q2NVfXSgVQlSZpx+g2U0wZZhCRp5usrUKrqyiRPAw5omq6qqrsGV5Ykaabpaw4lyZ8DVwGvAP4c+F6SYwdZmCRpZun3ktcpwAEbzkqSzAUuBT47qMIkSTNLv+/y2mrMJa57NuGxkqQtQL9nKF9LcjHw6Wb9OOArgylJkjQT9Tsp//YkLwf+BAiwvKo+P9DKJEkzSt+f5VVVnwM+N8BaJEkz2EYDJcm3qurgJPfT+eyuxzcBVVVPHWh1kqQZY6OBUlUHN/fbT085kqSZqt+/Qzm/n7ZNlWRJkuuTrEmyrMf2OUkubLZ/L8nCrm0nNe3XJzlyqrVIkqam37f+7te9kmQ28LypHDjJLOBDwFHAYuCVSRaP6fZa4N6qWgScA5zdPHYxcHxT1xLgw83+JElDstFAac4C7geeneSXze1+4E7gi1M89oHAmqq6saoeAS4Ajh7T52jgvGb5s8ChSdK0X1BVD1fVTcCaZn+SpCHZaKBU1ZnN/Mn7quqpzW37qtqlqk6a4rHnAbd2ra9t2nr2qapHgfuAXfp8LABJTkwymmR03bp1UyxZkjSefv8O5aQkOwH7Alt3tX9jCsdOr0P12aefx3Yaq5YDywFGRkb8gjBJGpB+vwL4dcCbgPl0vgb4PwLfAV40hWOvBfbqWp8P3D5On7XNvM0OwPo+HytJmkb9Tsq/ic5H199SVS8EngNM9frR1cC+SfZJ8mQ6k+wrxvRZASxtlo8FLm++hngFcHzzLrB96Jw5XTXFeiRJU9DvX8o/VFUPJSHJnKr6aZJnTuXAVfVokjcCFwOzgI9X1aokZwCjVbUC+BhwfpI1dM5Mjm8euyrJRcCPgUeBN1TVY1OpR5I0Nf0GytokOwJfAC5Jci8tXGKqqq8w5kMmq+pvu5YfovMdLL0e+x7gPVOtQZLUjn4n5f+sWTwtydfpzGV8bWBVSZJmnAkDJclWwA+r6lnQ+TrggVclSZpxJpyUr6rfAj9IsmAa6pEkzVD9zqHsAaxKchXwqw2NVfXSgVQlSZpx+g2U0wdahSRpxut3Uv7KJHsD+1bVpUm2pfNWX0mSgP4/vv71dD6c8f80TfPovIVYkiSg/7+UfwOd75P/JUBVrQZ2G1RRkqSZp99Aebj5iHng8e9D8YMWJUmP6zdQrkxyMrBNksOBzwBfGlxZkqSZpt9AWUbnwyB/BPwl8JWqOmVgVUmSZpx+3zb8V1X198BHNzQkeVPTJklS32coS3u0ndBiHZKkGW6jZyhJXgm8CtgnSfd3lWwP3DPIwiRJM8tEl7y+DdwB7Aq8v6v9fuCHgypKkjTzbDRQquoW4BbgoOkpR5I0U010yetbVXVwkvt54t+dBKiqeupAq5MkzRgTnaEc3NxvPz3lSJJmqn7f5SVJ0kYZKJKkVhgokqRWGCiSpFYYKJKkVhgokqRWGCiSpFYYKJKkVhgokqRWGCiSpFYYKJKkVhgokqRWDCVQkuyc5JIkq5v7ncbpt7TpszrJ0qZt2yRfTvLTJKuSnDW91UuSehnWGcoy4LKq2he4rFl/giQ7A6cCzwcOBE7tCp6/q6r/ADwH+JMkR01P2ZKk8QwrUI4GzmuWzwOO6dHnSOCSqlpfVfcClwBLqurBqvo6QFU9AlwDzJ+GmiVJGzGsQHlaVd0B0Nzv1qPPPODWrvW1TdvjkuwIvITOWY4kaYgm+k75SUtyKbB7j02n9LuLHm2Pf2tkktnAp4F/qKobN1LHicCJAAsWLOjz0JKkTTWwQKmqw8bbluTOJHtU1R1J9gDu6tFtLXBI1/p84Iqu9eXA6qr6wAR1LG/6MjIyUhvrK0mavGFd8loBLG2WlwJf7NHnYuCIJDs1k/FHNG0keTewA/DmaahVktSHYQXKWcDhSVYDhzfrJBlJci5AVa0H3gVc3dzOqKr1SebTuWy2GLgmybVJXjeMQUiSfidVW85VoJGRkRodHR12GZI0oyRZWVUjE/XzL+UlSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0wUCRJrTBQJEmtMFAkSa0YSqAk2TnJJUlWN/c7jdNvadNndZKlPbavSHLd4CuWJE1kWGcoy4DLqmpf4LJm/QmS7AycCjwfOBA4tTt4krwMeGB6ypUkTWRYgXI0cF6zfB5wTI8+RwKXVNX6qroXuARYApBkO+CtwLunoVZJUh+GFShPq6o7AJr73Xr0mQfc2rW+tmkDeBfwfuDBiQ6U5MQko0lG161bN7WqJUnjmj2oHSe5FNi9x6ZT+t1Fj7ZKsj+wqKrekmThRDupquXAcoCRkZHq89iSpE00sECpqsPG25bkziR7VNUdSfYA7urRbS1wSNf6fOAK4CDgeUluplP/bkmuqKpDkCQNzbAuea0ANrxraynwxR59LgaOSLJTMxl/BHBxVf3vqtqzqhYCBwP/bphI0vANK1DOAg5Psho4vFknyUiScwGqaj2duZKrm9sZTZskaTOUqi1nWmFkZKRGR0eHXYYkzShJVlbVyET9/Et5SVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSKwwUSVIrDBRJUisMFElSK1JVw65h2iRZB9wy7Do20a7A3cMuYpo55i2DY5459q6quRN12qICZSZKMlpVI8OuYzo55i2DY/794yUvSVIrDBRJUisMlM3f8mEXMASOecvgmH/POIciSWqFZyiSpFYYKJKkVhgom4EkOye5JMnq5n6ncfotbfqsTrK0x/YVSa4bfMVTN5UxJ9k2yZeT/DTJqiRnTW/1mybJkiTXJ1mTZFmP7XOSXNhs/16ShV3bTmrar09y5HTWPRWTHXOSw5OsTPKj5v5F0137ZEzlZ9xsX5DkgSRvm66aB6KqvA35BrwXWNYsLwPO7tFnZ+DG5n6nZnmnru0vAz4FXDfs8Qx6zMC2wAubPk8GvgkcNewxjTPOWcANwNObWn8ALB7T578DH2mWjwcubJYXN/3nAPs0+5k17DENeMzPAfZslp8F3Dbs8QxyvF3bPwd8BnjbsMczlZtnKJuHo4HzmuXzgGN69DkSuKSq1lfVvcAlwBKAJNsBbwXePQ21tmXSY66qB6vq6wBV9QhwDTB/GmqejAOBNVV1Y1PrBXTG3q37ufgscGiSNO0XVNXDVXUTsKbZ3+Zu0mOuqu9X1e1N+ypg6yRzpqXqyZvKz5gkx9D5ZWnVNNU7MAbK5uFpVXUHQHO/W48+84Bbu9bXNm0A7wLeDzw4yCJbNtUxA5BkR+AlwGUDqnOqJhxDd5+qehS4D9ilz8dujqYy5m4vB75fVQ8PqM62THq8SZ4CvAM4fRrqHLjZwy5gS5HkUmD3HptO6XcXPdoqyf7Aoqp6y9jrssM2qDF37X828GngH6rqxk2vcFpsdAwT9OnnsZujqYy5szHZDzgbOKLFugZlKuM9HTinqh5oTlhmNANlmlTVYeNtS3Jnkj2q6o4kewB39ei2Fjika30+cAVwEPC8JDfT+XnuluSKqjqEIRvgmDdYDqyuqg+0UO6grAX26lqfD9w+Tp+1TUjuAKzv87Gbo6mMmSTzgc8Dr66qGwZf7pRNZbzPB45N8l5gR+C3SR6qqg8OvuwBGPYkjrcCeB9PnKB+b48+OwM30ZmU3qlZ3nlMn4XMnEn5KY2ZznzR54Cthj2WCcY5m8718X343YTtfmP6vIEnTthe1CzvxxMn5W9kZkzKT2XMOzb9Xz7scUzHeMf0OY0ZPik/9AK8FXSuHV8GrG7uN7xojgDndvV7DZ2J2TXAf+2xn5kUKJMeM53fAAv4CXBtc3vdsMe0kbG+GPh3Ou8EOqVpOwN4abO8NZ13+KwBrgKe3vXYU5rHXc9m+k62NscMvBP4VdfP9Vpgt2GPZ5A/4659zPhA8aNXJEmt8F1ekqRWGCiSpFYYKJKkVhgokqRWGCiSpFYYKNKAJdk9yQVJbkjy4yRfSfIHm7iPkwdVn9QW3zYsDVDzAYDfBs6rqo80bfsD21fVNzdhPw9U1XYDKlNqhWco0mC9EPjNhjABqKprgW8leV+S65rv/jgOIMkeSb6R5Npm2wua73vZpmn75JDGIU3Iz/KSButZwMoe7S8D9gf+GNgVuDrJN4BXARdX1XuSzAK2rapvJnljVe0/bVVLk2CgSMNxMPDpqnoMuDPJlcABwNXAx5M8CfhCczYjzQhe8pIGaxXwvB7tPT+rvKq+AfwpcBtwfpJXD7A2qVUGijRYlwNzkrx+Q0OSA4B7geOSzEoyl06IXJVkb+Cuqvoo8DHguc3DftOctUibLS95SQNUVZXkz4APJFkGPATcDLwZ2I7OR50X8DdV9fMkS4G3J/kN8ACw4QxlOfDDJNdU1V9M9zikfvi2YUlSK7zkJUlqhYEiSWqFgSJJaoWBIklqhYEiSWqFgSJJaoWBIklqxf8DMnn/1+XlIcEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"stochastic\",\"train\")\n",
    "plt.title(\"stochastic gradient descent\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![4.2 derivation](42updated.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The difference in converge between batch and stochastic gradient descent is that stochastic converges faster (and more consistently between different runs)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 stochastic gradient descent with decaying learning rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/cs189/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEWCAYAAADLkvgyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecFdX5x/HPF5belrIgLFVZRECkrL2CDRMVNRI1JKLxF6MRe2KJSSwpmhijUaOJ0ShWJCZGYgmiFCvK0rssSFlAWHrvz++POavXdctl2WXu7j7v1+u+9s6ZmTPP3J07z52ZM3NkZjjnnHNxqhF3AM4555wnI+ecc7HzZOSccy52noycc87FzpORc8652Hkycs45FztPRilEkknqfACW85akIRW9nH0l6RlJvwnvT5Q0L8ZYDsj/Ik6SZkk6pYTx4yT9337UH9tnKKm9pM2SalZA3V9upwdaqn53y4Mno/0k6S5Jz8cdR3GKis/MzjKzYXHFlAwze9/MDi2PuiQtknRaedQVN0mXSfqgPOoys+5mNi7Um9Lb8b4ysyVm1tDM9sQdS3lKpe/u/v5YKcyTkasQktLijsG5yiSVvjOxxGJm/kriBdwKLAM2AfOAU4EBwE5gF7AZmBambQOMBNYCucCPEuqpCfwcWBDqmgS0C+MMuAqYD6wD/gIojDsEGAOsAVYDLwDpZYxvHPB/CfP+CJgT5p0N9CnmMzgj1L0BeAwYX1APcBnwIfBgWO/fJBFzb2ByWO7LwHDgN2HcKUBewrRtgH8B+cDnwHUJ4+4CRgDPhrpmAdlh3HPAXmBb+AxuKWbdfgasAJYDPwz/i85hXB3gj8ASYCXwV6BeGNcCeB1YH9b7faBGGNcO+HeIeQ3waMLyfhg+83XAKKBDwrgitwPgMGA7sCesy/oi1qMfMCNh+B3g04ThD4DzwvtFwGmlbCe/Dv/XTcDbQIsSviNl+gzD+IHAVGAj0XdjQCi/nK+2zYXAjxPmmQmckzBci2g76wV0DMtPS2ZdgEuBxeH/9MuCz6aY9XyGsJ2G4bND7OuBj4CeCeNu46vv+mzg/IRxl/HN78xl4X/0x/C//xw4K2GecXz9O1fStJ2A98Ky3yHajp4vZp1OAfKI9iNfEH1vmhJt2/mh/teBtmH63xJth9uJtplHQ3lXYHRYn3nAd5Pex8a5g68sL+BQYCnQJgx3BA4J7+8q/A8m2kk/BtQNX4x84NSEL+yMUKeAI4DmYZyFf3g60D7MV/Cl7AycTvSlzggb2UNljC9xgx5ElMSODPF0JmHHmDBPC6IdxQVAGnA90c4r8YuxG7g2jK9XSsy1ib78NxLtRC4M9X0jGREdwU8CfhXmO5hox3RmwjpuB75FlOzvBSYkxL6IYnYsYfwAoh1kD6AB8CJf35E+RPTjohnQCPgvcG8Ydy/RjrVWeJ0YPseawDSiHU2DsC2cEOY5j+hHymHhs/oF8FFCPCVtB5cBH5SwLnWJEm+LUPcXRMmhUfifbOOr7e3Lz4Xit5MFQJcw7zjgvgr4DI8i+oFzevhfZwJdw7hvE/2oEXAysJXwYwm4BXg5IYaBhERM0cmoyHUBuhHtUE8g2r7+SLQtlpqMgD7AKuDo8D8fEj7XOgnfrzZhvS4CtgCtS/jOXBaW/aNQ39Xh/6eE9Uj8zpU07cdhXWqHddtY+H+csE6nhFh+T/R9rQc0B74D1A//s38C/ylqPxKGGxDthy4P69OH6MdB96T2s3Hv6CvDi2inuoroV2StQuPuSvwHE/0a3gM0Sii7F3gmvJ8HDCxmOUbYYYXhEcBtxUx7HjBlX+MrYoMeBVyfxGdwKfBxwrDChpf4xVhSSh2JMZ+U+MUJZR9RdDI6unDdwO3A0wnr+E7CuG7AtoThRZScjP5Bwk6WaIdl4XMV0Q7kkITxxwKfh/f3AK8RdrqFpskn7AwLjXsLuCJhuAbRTrZDadsBpSSjMM37RD8ajiE6AhhBlCz6AdOL+lxK2E5+kTD8E+B/FfAZ/g14MMnv4n8KtleinfwmoHEYfoVw5EvRyajIdSH6kfNSwrj6REeKySSjx4FfFxo/Dzi5mHmnEr7/FPGdCWW5hWIx4KAivrvFTkv0I2Y3UD9h/POF/8cJ404J61y3hM++F7Cu0PaRmIwuAt4vNM/fgDuT+d/6NaMkmFkucAPRF3aVpOGS2hQzeRtgrZltSihbTPRrD6JktaCExX2R8H4r0BBAUsuw3GWSNhJtWC3KEF9hpcVToA1R8iEs04gO6xMtTRwoKeZQ37JQT4HFxSy7A9BG0vqCF9GpzlYJ0xT+3Oruw3nvr61boTgyiL7kkxKW/b9QDnA/0VHO25IWSrotlLcDFpvZ7mLW588J9a0l2mFnJkxT5HaQpPFEO5eTwvtxREcVJ4fhfZFsHPvzGRa7DUo6S9IESWvDfN/iq+1+OdFpru9ISgfOIjoVvK/rUnjb3kp0ui4ZHYCbC22b7UKdSLpU0tSEcT346jsAhb4zheMMsUDxn3tx0xbsh7YmTFvUshLlm9n2ggFJ9SX9TdLi8P19D0gvoYViB+DoQp/FYKLkWCpPRkkysxfN7ASiD9yIDmcJ7xMtB5pJapRQ1p7oVBhEG8QhZQjh3rCsnmbWGPg+0Q5sX+MrLNl4VgBtCwYkKXG4mGWVFPMKIDPUU6B9CTF+bmbpCa9GZvatJOIuKq7CVhDtQIqKYzXRqa3uCctuYmYNAcxsk5ndbGYHA+cAN0k6NcTcvpiEuJTo2kfi+tQzs4/KYV3gm8loPKUno2TqLUmZP0OK2QYl1SG6TvhHoJWZpQNvkrDdA8OItqtBREfuywrXk2Tsidt2wSmqZCwFflvof1nfzF6S1AH4OzCU6NRoOtF1rsT49/dzL84Kov1Q/YSydsVNXEwsNxNdAjg6fH9PCuUqZvqlwPhCn0VDM7s6mYA9GSVB0qGS+ocvx3aiL1ZBk9GVQEdJNQDMbCnR6aZ7JdWV1BO4gq9+sT0J/FpSliI9JSWz4TciXLSWlEl07Wmf4yvCk8BPJfUN8XQOX6LC3gAOl3Re2MFeQ+m/eIqNmeh89m7gOklpki4gunZQlE+BjZJulVRPUk1JPSQdWcryC6wkus5UnBHAZZK6hS/vnQUjzGwv0Q7lQUktASRlSjozvD87fGYiOie/J7w+Jdoh3CepQdgWjg/V/hW4XVL3UEcTSYP2YV3aSqpdwjQfEe1EjiJqvDCL8KuV6NdtcfWWtJ2UpsyfIfAUcLmkUyXVCOO6El3rqEN0unO3pLOIGtEk+g/RtYnriRqwlMUrwDmSjguf6918PWGU5O/AVZKODt+fBpK+HX6MNiDaYeeHdb6c6MiowpnZYiAHuEtSbUnHEv1Y2heNiPYl6yU1I+F/GhT+Xr0OdJH0A0m1wutISYclszBPRsmpA9xH9AvvC6Al0WkiiC7qAayRNDm8v4TonPVy4FWic6ajw7g/EX1x3ybaeT1FdLGwNHcTfek2ECWGf+9HfF8ys38StYx5kej8+3+ILjIXnm410a/PPxCdwuhGtLHvKEvMZraT6LrGZUQtdS4qtE6Jy95D9EXqRdRiaDVREm1SwrIT3Qv8Ipw6+GkR9b9FdIF9DNEptzGFJrk1lE8IpyveIdrZA2SF4c1ECfYxMxuXEHNnohZkeWEdMbNXiY5ch4f6ZhKdYkrGGKLWgl9IWl3UBGa2haiV4qzwORNiW2xmq4qpt8TtpDT78xma2adEF70fJNpWxhNdP9sEXEf0fVkHfI+oEUTicrcRHT11opjtJ4nYZxE1IhhO9ANiE9E12JK27YJ5c4gaEDwaYswl2qYxs9nAA0Sf/UrgcKLTigfKYKJrc2uIWuq9TBLrlOAhon3TamAC0anVRH8GLpS0TtLD4f91BnAx0b7vC75qEFGqglYXzu2T8As6DxhsZmPjjsdVX5J+BXQxs++XU30NiZppZ5nZ5+VRZyqQ9DIw18wKH+GkBD8yckmTdKak9HA68OdEpzImxByWq8bC6aMrgCf2s55zwgX7BkTXqGYQtTastMIpskPCqc8BRE3f/xN3XMXxZOT2xbFErZ5WE52COi+cJnHugJP0I6KL5m+ZWXHXwpI1kOjU0nKiU68XW+U/bXQQUUvKzcDDwNVmNiXWiErgp+mcc87Fzo+MnHPOxS5lHsyX6lq0aGEdO3aMOwznnKtUJk2atNrMMkqbzpNRkjp27EhOTk7cYTjnXKUiqbgnq3yNn6ZzzjkXO09GzjnnYufJyDnnXOw8GTnnnIudJyPnnHOx82TknHMudp6MnHPOxc6TUQXbuH0Xd42cxYZtu+IOxTnnUpYnowq2YNVmnp+wmFtemYY/B9A554rmyaiC9W7flFsGHMqoWSt5+sNFcYfjnHMpyZPRAfCjEw/mtMNa8rs35zBlybq4w3HOuZQTSzKSNEjSLEl7JWUXMb69pM2JXURLGiBpnqRcSbcllHeS9Imk+ZJeDn3YI6lOGM4N4zsmzHN7KJ8n6cyKXVuQxAODetGqcV2GvjiF9Vt3lj6Tc85VI3EdGc0ELgCK6xDrQeCtggFJNYG/AGcB3YBLJHULo38PPGhmWUR90F8Ryq8A1plZ51Df70Nd3Yj6aO8ODAAeC/VXqCb1a/GXwX1YtWk7N4+Yxt69fv3IOecKxJKMzGyOmc0rapyk84CFwKyE4qOAXDNbaGY7geHAQEkC+gOvhOmGAeeF9wPDMGH8qWH6gcBwM9sR+rfPDfVXuF7t0vn5tw7j3bmr+Pv7Cw/EIp1zrlJIqWtGof/5W4G7C43KJOpeuEBeKGsOrDez3YXKvzZPGL8hTF9cXUXFc6WkHEk5+fn5ZV2tr7nsuI6c1eMg/jBqHjmL1pZLnc45V9lVWDKS9I6kmUW8BpYw291Ep9w2F66uiGmthPKyzvP1QrMnzCzbzLIzMkrtGyopkvj9hT3JTK/H0BensGbzjnKp1znnKrMK61zPzE4rw2xHAxdK+gOQDuyVtB2YBLRLmK4tsBxYDaRLSgtHPwXlEB3xtAPyJKUBTYC1CeWF6zpgGtetxWOD+3DBYx9x44hpPHPZkdSoUVSOdM656iGlTtOZ2Ylm1tHMOgIPAb8zs0eBiUBWaDlXm6gBwkiL7iIdC1wYqhgCvBbejwzDhPFjwvQjgYtDa7tOQBbw6QFYva/pkdmEX53Tjfc+y+fx8QsO9OKdcy6lxNW0+3xJecCxwBuSRpU0fTjqGQqMAuYAI8ysoIHDrcBNknKJrgk9FcqfApqH8puA20Jds4ARwGzgf8A1ZranPNcvWYOPbs85R7Thgbfn8fGCNXGE4JxzKUH+iJrkZGdnW05OTrnXu3nHbs595AM27djNm9edSEajOuW+DOeci4ukSWb2jftJC0up03TVUcM6afxlcB82btvF9cOnsMfvP3LOVUOejFLAYa0bc8/A7ny0YA0Pvzs/7nCcc+6A82SUIr6b3Y4L+mTy8Jj5fDB/ddzhOOfcAeXJKEVI4jfn9aBzRkOuHz6FlRu3xx2Sc84dMJ6MUkj92mk8NrgPW3fu4dqXprB7z964Q3LOuQPCk1GKyWrViN+e34NPP1/Ln0Z/Fnc4zjl3QHgySkEX9GnLRdnteGzcAsbOWxV3OM45V+E8GaWouwd2p+tBjbjp5aksX78t7nCcc65CeTJKUXVr1eSxwX3YuXsvQ1+czC6/fuScq8I8GaWwgzMacu93ejJ5yXruH1Vk90/OOVcleDJKcece0YbvH9OeJ95byOjZK+MOxznnKoQno0rgF9/uRvc2jbl5xFSWrt0adzjOOVfuPBlVAgXXj8xg6IuT2bnbrx8556oWT0aVRIfmDfjDhT2ZlreB3705J+5wnHOuXHkyqkTOOrw1lx3XkWc+WsRbM1bEHY5zzpWbuDrXGyRplqS9krITyjtK2iZpanj9NWFcX0kzJOVKeliSQnkzSaMlzQ9/m4ZyhelyJU2X1CehriFh+vmShlCJ/Pxbh3FEu3RueWU6i9dsiTsc55wrF3EdGc0ELgDeK2LcAjPrFV5XJZQ/DlxJ1E14FjAglN8GvGtmWcC7YRjgrIRprwzzI6kZcCdwNHAUcGdBAqsMaqfV4NFLeiPBT16YzPZdsXRS65xz5SqWZGRmc8ws6RtnJLUGGpvZxxZ1TfsscF4YPRAYFt4PK1T+rEUmAOmhnjOB0Wa21szWAaP5KrFVCu2a1eeB7/Zi1vKN/OaN2XGH45xz+y0Vrxl1kjRF0nhJJ4ayTCAvYZq8UAbQysxWAIS/LRPmWVrEPMWVf4OkKyXlSMrJz8/fn3Uqd6d3a8WVJx3M8xOW8NrUZXGH45xz+6XCkpGkdyTNLOI1sITZVgDtzaw3cBPwoqTGgIqYtrT+uYubJ+m6zOwJM8s2s+yMjIxSFnfg/ezMQ+nboSk///cMFuRvjjsc55wrswpLRmZ2mpn1KOL1Wgnz7DCzNeH9JGAB0IXo6KVtwqRtgeXh/cpw+q3gdF7BY67zgHZFzFNceaVTq2YNHrmkN7XTanCNXz9yzlViKXWaTlKGpJrh/cFEjQ8WhtNvmyQdE1rRXQoUJLWRQEGLuCGFyi8NreqOATaEekYBZ0hqGhounBHKKqU26fX400W9mPvFJu58bVbc4TjnXJnE1bT7fEl5wLHAG5IKksFJwHRJ04BXgKvMbG0YdzXwJJBLdMT0Vii/Dzhd0nzg9DAM8CawMEz/d+AnAKG+XwMTw+uehGVUSv0Obck1/Q7h5Zyl/GtSXukzOOdcilHUOM2VJjs723JycuIOo1i79+xl8JOfMD1vAyOHHk9Wq0Zxh+Scc0iaZGbZpU2XUqfpXNml1azBw5f0pn7tmvzkhcls3bk77pCccy5pnoyqkFaN6/Lni3uTm7+ZX7w6Ez/qdc5VFp6MqpgTslpwXf8s/j1lGSNylpY+g3POpQBPRlXQdadmcXzn5vzqtVnMWbEx7nCcc65UnoyqoJo1xEMX9aZxvVpc88JkNu/w60fOudTmyaiKymhUh0cu6c2iNVu4/d8z/PqRcy6leTKqwo45uDk3n3Eo/522nBc+WRJ3OM45VyxPRlXc1ScfwsldMrjnv7OZuWxD3OE451yRPBlVcTVqiAcv6kWzBrX5yQuT2bh9V9whOefcN3gyqgaaNajNo9/rzbL127j1lel+/cg5l3I8GVUT2R2bceuAQ3lr5hc889GiuMNxzrmv8WRUjfzoxIM57bCW/O7NOUxduj7ucJxz7kuejKoRSfxx0BG0bFSXa16YzPqtO+MOyTnnAE9G1U56/ej60apN2/npP/36kXMuNXgyqoZ6t2/K7WcdxjtzVvLk+5/HHY5zzsXWud4gSbMk7ZWUXWhcT0kfh/EzJNUN5X3DcK6kh0OPr0hqJmm0pPnhb9NQrjBdrqTpkvokLGNImH6+pCFUQ5cf35EB3Q/ivv/NZdLiSt23oHOuCojryGgmcAHwXmKhpDTgeaIeXrsDpwAFN8Y8DlxJ1BV5FjAglN8GvGtmWcC7YRjgrIRprwzzI6kZcCdwNHAUcGdBAqtOJPGHQT3JTK/H0BensHaLXz9yzsUnlmRkZnPMbF4Ro84AppvZtDDdGjPbI6k10NjMPrboIsezwHlhnoHAsPB+WKHyZy0yAUgP9ZwJjDaztWa2DhjNV4mtWmlctxaPDe7Dms07ufHlqezd69ePnHPxSLVrRl0AkzRK0mRJt4TyTCAvYbq8UAbQysxWAIS/LRPmWVrEPMWVf4OkKyXlSMrJz8/fj9VKXT0ym/DLc7ox/rN8Hh+/IO5wnHPVVFpFVSzpHeCgIkbdYWavlRDPCcCRwFbgXUmTgKI65SntZ7yKmae48m8Wmj0BPAGQnZ1dZQ8bvn90ez5ZuIYH3p5H3w5NOebg5nGH5JyrZirsyMjMTjOzHkW8iktEEB2ljDez1Wa2FXgT6BPK2yZM1xZYHt6vDKffCH9XJdTVroh5iiuvtiRx33d60rF5A657aQr5m3bEHZJzrppJtdN0o4CekuqHxgwnA7PD6bdNko4JreguBQqS2kigoEXckELll4ZWdccAG0I9o4AzJDUNDRfOCGXVWsM6afxlcB82bNvFjS9PZY9fP3LOHUBxNe0+X1IecCzwhqRRAKFBwZ+AicBUYLKZvRFmuxp4EsgFFgBvhfL7gNMlzQdOD8MQHVUtDNP/HfhJWMZa4NdhGROBe0JZtXdY68bcfW53PshdzSNj5scdjnOuGpHfgZ+c7Oxsy8nJiTuMCmdm3DxiGq9OXcbzVxzN8Z1bxB2Sc64SkzTJzLJLmy7VTtO5mEniN+f34JCMhlw/fAqrNm6POyTnXDXgych9Q/3aaTw+uA9bduzh2pemsHvP3rhDcs5VcZ6MXJGyWjXiN+f14JPP1/LQO379yDlXsTwZuWJ9p29bvpvdlkfH5jJu3qrSZ3DOuTLyZORKdPe5Peh6UCNufHkqKzZsizsc51wV5cnIlahe7Zr8ZXAfdu7ey7UvTmGXXz9yzlUAT0auVIdkNOR3FxxOzuJ1/HFUUc+3dc65/ePJyCVlYK9MBh/dnr+9t5D/zVwRdzjOuSrGk5FL2q/O6cYR7dK5ecQ0cldtjjsc51wV4snIJa1OWk0eH9yHurVq8uPncti8Y3fcITnnqghPRm6ftEmvxyPf682iNVv56Yhp+OOknHPlwZOR22fHHdKC28/qyv9mfcFfxy+MOxznXBXgyciVyRUndOLsnq25f9Rc3p9fNXvBdc4dOJ6MXJlI4g8X9iSrZSOue2kKS9dujTsk51wl5snIlVn92mn89Qd92b3XuPqFSWzftSfukJxzlVRcnesNkjRL0l5J2QnlgyVNTXjtldQrjOsraYakXEkPhx5fkdRM0mhJ88PfpqFcYbpcSdMl9UlYzpAw/XxJQwrH55LXqUUDHrqoFzOXbeSOV2d6gwbnXJnEdWQ0E7gAeC+x0MxeMLNeZtYL+AGwyMymhtGPA1cCWeE1IJTfBrxrZlnAu2EY4KyEaa8M8yOpGXAncDRwFHBnQQJzZXPqYa247tQs/jU5j+c/WRJ3OM65SiiWZGRmc8ystOfKXAK8BCCpNdDYzD626Kf3s8B5YbqBwLDwflih8mctMgFID/WcCYw2s7Whm/PRfJXYXBndcGoW/Q7N4J7/zmLSYu/F3Tm3b1L5mtFFhGQEZAJ5CePyQhlAKzNbARD+tkyYZ2kR8xRX7vZDjRrioYt607pJPa5+fjKrNnkPsc655FVYMpL0jqSZRbwGJjHv0cBWM5tZUFTEZKVdnChunqTrknSlpBxJOfn53ny5NE3q1+JvP+jLxu27uOaFyf6Eb+dc0iosGZnZaWbWo4jXa0nMfjFfHRVBdPTSNmG4LbA8vF8ZTr8VnM5blTBPuyLmKa68qHV4wsyyzSw7IyMjibDdYa0b8/vv9GTionX89o05cYfjnKskkk5Gko6T9D1Jlxa8KiIgSTWAQcDwgrJw+m2TpGNCK7pLgYKkNhIoaBE3pFD5paFV3THAhlDPKOAMSU1Dw4UzQpkrJwN7ZfLD4zvxzEeL+M+UZXGH45yrBNKSmUjSc8AhwFSg4GaSgoYE+0zS+cAjQAbwhqSpZnZmGH0SkGdmhZ8zczXwDFAPeCu8AO4DRki6AlhClMgA3gS+BeQCW4HLAcxsraRfAxPDdPeYmV9xL2e3f6srM5dv4LZ/T6dLq0Z0a9M47pCccylMydwXImkO0M2q8U0k2dnZlpOTE3cYlUr+ph2c/cj71E6rwX+HnkB6/dpxh+ScO8AkTTKz7NKmS/Y03UzgoP0LyVU3GY3q8Pj3+/LFhu1cP3wqe/ZW298yzrlSJJuMWgCzJY2SNLLgVZGBuaqhT/um3HVud8Z/ls+f3/ks7nCccykqqWtGwF0VGYSr2r53VHumLV3Pw2NyObxtOqd3axV3SM65FJPUkZGZjQfmAo3Ca04oc65UkrhnYA96tm3CTS9PZWG+d1nunPu6pJKRpO8CnxK1VPsu8ImkCysyMFe11K1Vk8cG9yGtpvjxc5PY4l2WO+cSJHvN6A7gSDMbYmaXEj1g9JcVF5arito2rc8jl/RhQf5mbnlluj/h2zn3pWSTUQ0zW5UwvGYf5nXuSydkteCWAV15Y8YK/v6+d1nunIsk24Dhf5JG8dUjei4iuqnUuX3245MOZtrS9dz31lx6tGnCcZ1bxB2Scy5myTZg+BnwBNATOAJ4wsxurcjAXNUlifsHHcHBGQ0Z+tIUlq/fFndIzrmYJX2qzcz+ZWY3mdmNZvZqRQblqr6GddL42w/6snP3Xq5+3rssd666KzEZSfog/N0kaWPCa5OkjQcmRFdVHZLRkAe+ewTT8jZw18hZcYfjnItRicnIzE4IfxuZWeOEVyMz8ydfuv12ZveDGNqvM8MnLuWlT73Lcueqq2TvM3oumTLnyuLG07twUpcM7nxtFlOWrIs7HOdcDJK9ZtQ9cUBSGtC3/MNx1VHNGuLhi3vRqkkdrn5+MvmbdsQdknPuACvtmtHtkjYBPROvFwEr+aoTO+f2W3r92vz1+31Zt3UnQ1+czG7vsty5aqW0a0b3mlkj4P5C14uam9ntByhGV010b9OE+75zOJ98vpb73pobdzjOuQMo2fuMbg/ddB8l6aSCV1kXKmmQpFmS9krKTiivJWmYpBmS5ki6PWHcAEnzJOVKui2hvJOkTyTNl/SypNqhvE4Yzg3jOybMc3sonyepoIdZlwLO792WIcd24MkPPmfktOVxh+OcO0CSbcDwf8B7wCjg7vD3rv1Y7kzgglBnokFAHTM7nOia1I8ldZRUE/gLcBbQDbhEUrcwz++BB80sC1gHXBHKrwDWmVln4MEwHWG+i4mugw0AHgv1uxRxx7e7kd2hKbe+Mp25X/gdBM5VB8k2YLgeOBJYbGb9gN5AflkXamZzzGxeUaOABqGBRD1gJ7CR6MGsuWa20Mx2AsOBgZIE9AdeCfMPA84L7weGYcL4U8P0A4HhZrbDzD4HckP9LkXUTqvBY4P70LBuGj9+bhIbtu2KOyTnXAVLNhltN7PtEJ3+MrO5wKEVEM8rwBZgBbAE+KOZrQUygaUJ0+WFsubAejPbXaicxHnC+A1h+uLq+gZJV0rKkZSTn1/m3OvKoGXjujw+uA/L1m3jppeUnuQVAAAdZ0lEQVSnste7LHeuSks2GeVJSgf+A4yW9BpQ4gl9Se9ImlnEa2AJsx0F7AHaAJ2AmyUdDKiIaa2Ecso4z9cLzZ4ws2wzy87IyCghbFcRsjs241fndOPduat4ZExu3OE45ypQUk/tNrPzw9u7JI0FmgD/K2We08oQz/eA/5nZLmCVpA+BbKIjmXYJ07UlSoargXRJaeHop6AcoiOedkSJNC3EvDahvHBdLgX94JgOTF26nofe/YzD2zamf1fvsty5qqjUIyNJNSTNLBg2s/FmNjJcuylvS4D+ijQAjiHq7nwikBVaztUmaoAw0qLe2cYCBb3ODuGr+59GhmHC+DFh+pHAxaG1XScgi6gXW5eCJPG78w/nsIMac8PwqSxavSXukJxzFaDUZGRme4FpktqX10IlnS8pDzgWeCP0lQRRi7mGRK3tJgJPm9n0cNQzlKgV3xxghJkVPFnzVuAmSblE14SeCuVPAc1D+U3AbWF9ZgEjgNlER3fXmJk/MjqF1a1Vk7/9oC81aoirnp/E1p3eZblzVY2S6fpZ0hii1nSfEjUwAMDMzq240FJLdna25eTkxB1GtfbeZ/kMefpTzunZhj9f3IuocaRzLpVJmmRm2aVNl2xPr3fvZzzO7beTumTw0zMO5f5R8ziiXTpXnNAp7pCcc+Uk2QYM4yV1ALLM7B1J9QG/UdQdcD855RCm563nd2/OoXubxhxzcPO4Q3LOlYNkn8DwI6J7gP4WijKJmnk7d0BJ4o+DjqBD8/oMfXEyKzZ4l+XOVQXJ3md0DXA80dMQMLP5QMuKCsq5kjSqW4snftCXbTv3cPXzk9mx29ufOFfZJZuMdiQ25Q737fgt8S42nVs24o+DjmDq0vXc89/ZcYfjnNtPySaj8ZJ+DtSTdDrwT+C/FReWc6U76/DWXHXyIbzwyRJGTFxa+gzOuZSVbDK6jejBqDOAHwNvmtkdFRaVc0n66RldOKFzC37x2kym562POxznXBklm4yuNbO/m9kgM7vQzP4u6foKjcy5JKTVrMHDl/Qmo2EdrnpuEms2e5flzlVGySajIUWUXVaOcThXZs0aRF2Wr96yk+uGT/Euy52rhEpMRpIukfRfoJOkkQmvscCaAxOic6U7vG0TfnteDz7MXcP9bxfVVZZzLpWVdtPrR0R9C7UAHkgo3wRMr6ignCuLQdntmJa3nr+NX0jPzHS+3bN13CE555JUYjIys8XAYqIHmjqX8n51dndmLd/Iz16ZRpdWDclq1SjukJxzSSjtNN0H4e8mSRsTXpskbTwwITqXvNppNXh8cF/q1466LN+43bssd64yKDEZmdkJ4W8jM2uc8GpkZo0PTIjO7ZuDmtTlscF9WLJ2KzePmOZdljtXCSTbms65SuWoTs34+bcOY/TslTw2zrssdy7VxZKMJA2SNEvSXknZCeW1JT0taYakaZJOSRjXN5TnSnpYoTMbSc0kjZY0P/xtGsoVpsuVNF1Sn4S6hoTp50sqqtm6qwIuP74jA3u14YHRnzH+s/y4w6kwO3fvJZl+yZxLZcn2Z1TeZgIX8NVTwAv8CMDMDpfUEnhL0pGht9nHgSuBCcCbwADgLaKnQ7xrZvdJui0M3wqcRdSleBZwdJj/aEnNgDuBbKLn602SNNLM1lXkCrsDTxL3XnA4877YxHUvTeH1a0+gXbP6cYe1T8yM1Zt3snz9Npav38aygte6bSzfsI3l67ezdstO2jerT/+uLenftSVHH9yMOmnew4urXGJJRmY2Byiqp85uwLthmlWS1gPZkpYCjc3s4zDfs8B5RMloIHBKmH8YMI4oGQ0EnrXoJ+MESemSWodpR5vZ2lDXaKLE9lJFrKuLV/3aafztB30555EP+PFzk/jX1cdRr3bq7Kh37N7DivXbWb5+G3kh4RQkneXrt7Ns/TZ27v76Tbz1a9ckM70emU3r0bNtOhkN6zBj2QZe+nQJz3y0iPq1a3J85xb079qSfoe25KAmdWNaO+eSF9eRUXGmAQMlDQfaAX3D371AXsJ0eUR9KgG0MrMVAGa2IhxREcYvLWKe4sq/QdKVREdjtG/fvuxr5WLVoXkD/nxJb374zETueHUGD3z3iAPSZbmZsX7rrpBYtn3t77KQgPI3ffPxRS0b1SGzaT26tWnM6d1akZlejzbp9WiTXpe26fVpXC+tyPi379rDxwvWMGbuKsbMXcXo2SsB6Na6cZSYurakV7t0atbw7tpd6qmwZCTpHeCgIkbdYWavFTPbP4DDgByi+5s+AnYDRX17SjtJXtw8SddlZk8ATwBkZ2f7SflKrN+hLbnxtC78afRnHNEunSHHddzvOnft2csXG7YXSjTbv3y/fP02tu78el9LddJqfHlU0/XQlmQ2/XqiadWkTplPsdWtVZN+IencY8b8VZujxDRnFY+PX8CjY3Np1qA2J3fJoF/XlpyclUGT+rX2+3NwrjxUWDIys9PKMM9u4MaCYUkfAfOBdUDbhEnbAsvD+5WSWoejotbAqlCeR3RUVXiePL46rVdQPm5fY3WVz9B+nZmet55fvz6bbm0ac2THZiVOv3H7rq9Om637ZqJZuXE7hVuNN29Qm8ym9eic0ZCTsjLIbFqPzPS6tEmvR2Z6PZo1qH1Ajsok0aVVI7q0asRVJx/Chq27GD8/n7FzVzFu3ipenbKMmjVE3/ZN6ReuNXVp1fCAxOZcURRnKxxJ44CfmllOGK4fYtoS+k36pZmdFMZNBK4FPiFqwPCImb0p6X5gTUIDhmZmdoukbwNDgW8RNWB42MyOCg0YJgEFresmA30LriEVJzs723Jycsr3A3AH3Mbtuxj46Ids3rGbYZcfxbZdu8lbF12fKXwqbdP23V+bt1ZNRUcxTaKjmcKJpk16PerWSp3rUcXZs9eYunQ9Y8PpvNkrovvXM9Pr0a9rBv27tuS4Q1pUinVxqU/SJDPLLnW6OJKRpPOBR4AMYD0w1czOlNQRGEV0jWgZcEV4JBGhCfgzQD2ihgvXmplJag6MANoDS4BBZrY2NP1+lKhxwlbg8oSk90Pg5yGc35rZ06XF7Mmo6pj3xSbOf+zDb5xCa1Kv1pdJJTO9bsIptHq0Ta9Hi4Z1qFEFr7d8sWE7Y+dFienD3NVs3bmHOmk1OO6Q5l9ea2rbtHK1QnSpI6WTUWXkyahqmbV8A1OXrv8y0bROr0fDOqnWnufA27F7D58sXPtlI4gla7cC0KVVQ/p3bUX/ri3p0z6dtJp+v7xLjiejcubJyFU3ZsaC/C1fns6buGgtu/caTerV4qQuGfTvmsHJXVrSrEHtuEN1KcyTUTnzZOSqu43bd/HB/NWMCY0gVm/eiQS926V/eTqvW+vG3gjCfY0no3Lmyci5r+zda8xYtoExc1cxdt4qpudtAOCgxnXp1zWDfoe25PjOLWjgpz6rPU9G5cyTkXPFW7VpO+PmRU3H35+/ms07dlO7Zg2OPrjZl48p6tC8Qdxhuhh4MipnnoycS87O3XuZuChqBDF27ioWrt4CwMEZDeh/aEv6H9aSIzs2o5Y3gqgWPBmVM09GzpXNotVbvjyd98nCtezcs5dGddI4sUsL+h3aklMObUlGozpxh+kqiCejcubJyLn9t2XHbj7IXc3YkJxWboyezXdE2yb069qS0w5rRY/MJjFH6cqTJ6Ny5snIufJlZsxavjFqOj5vFVOXrscMzjmiDXee040WDf1oqSrwZFTOPBk5V7HWbN7BcxMW85exuTSok8Yvv92NC/pkelPxSi7ZZORXEJ1zKaF5wzrccFoX3rzuRA5u0YCb/zmNS//xKUvDUyBc1ebJyDmXUrJaNeKVq47j7nO7M3nxOs548D2e+uBz9hR+RLqrUjwZOedSTo0aYshxHXn7ppM5+uBm/Pr12Vzw+EfM/WJj3KG5CuLJyDmXsjLT6/H0ZUfy54t7sXTtVs5++AP+9PY8duzeU/rMrlLxZOScS2mSGNgrk3duOplzjmjDw2Ny+daf3ydnUYldkLlKxpORc65SaNagNg9e1ItnLj+S7bv2MuhvH/Or12ayecfu0md2KS+WZCTpfklzJU2X9Kqk9IRxt0vKlTRP0pkJ5QNCWW7o0bWgvJOkTyTNl/SypNqhvE4Yzg3jO5a2DOdc6jvl0Ja8feNJDDm2I89NWMzpfxrPmLkr4w7L7ae4joxGAz3MrCfwGXA7gKRuwMVAd6IeWh+TVFNSTeAvwFlAN+CSMC3A74EHzSwLWAdcEcqvANaZWWfgwTBdscuo4PV1zpWjBnXSuOvc7vzr6uNoVDeNHz6Tw3UvTWHN5h1xh+bKKJZkZGZvm1nBsfUEoG14PxAYbmY7zOxzIBc4KrxyzWyhme0EhgMDQ9fi/YFXwvzDgPMS6hoW3r8CnBqmL24ZzrlKpk/7prx+7YnccFoWb81cwWl/Gs+/J+fhN/NXPqlwzeiHwFvhfSawNGFcXigrrrw5sD4hsRWUf62uMH5DmL64ur5B0pWSciTl5Ofnl2nlnHMVq3ZajS9vlu3UogE3jZjGkKcnkrfOb5atTCosGUl6R9LMIl4DE6a5A9gNvFBQVERVVobystT1zUKzJ8ws28yyMzIyiprEOZcislo14p9XHcdd53QjZ9FaznjwPf7hN8tWGhXWDaOZnVbSeElDgLOBU+2rY+o8oF3CZG2B5eF9UeWrgXRJaeHoJ3H6grryJKUBTYC1pSzDOVeJ1awhLju+E6d3P4g7Xp3BPa/PZuS05fzhwp50adUo7vBcCeJqTTcAuBU418wSj6VHAheHlnCdgCzgU2AikBVaztUmaoAwMiSxscCFYf4hwGsJdQ0J7y8ExoTpi1uGc66KKLhZ9qGLerF4zRa+/fD7/Gn0Z36zbAqLq4P6R4E6wOjwRN4JZnaVmc2SNAKYTXT67hoz2wMgaSgwCqgJ/MPMZoW6bgWGS/oNMAV4KpQ/BTwnKZfoiOhigJKW4ZyrOiRxXu9MTsxqwa9fn83D787nzRkr+P13Dqdvh2Zxh+cK8S4kkuRdSDhXuY2dt4o7/j2DFRu3c+kxHfjZgK40rBPX7/Hqw7uQcM65BP0ObcnbN53MkGM78uyExZzxp/GMnbsq7rBc4MnIOVdtNAw3y75y1XE0qJPG5c9M5PrhfrNsKvBk5Jyrdvp2aMrr153ADadl8eaM6GbZV6f4zbJx8mTknKuW6qTV5IbTuvDGdSfSsUUDbnx5Gpf5zbKx8WTknKvWuoSeZe86pxsTw82yT3/oN8seaJ6MnHPVXsHNsm/feBJHdmzG3f+dzYV//YjPVm6KO7Rqw5ORc84FbZvW55nLo5tlF62ObpZ90G+WPSA8GTnnXIKCm2Xfuelkvn14a/787nzOfvgDJi1eF3doVZonI+ecK0LzhnV46OLePH3ZkWzZsZsL//oRd42c5T3LVhBPRs45V4J+Xb+6WXbYx4s488H3GDvPb5Ytb56MnHOuFF/dLHss9WrX5PKnJ3KD3yxbrjwZOedckvp2aMYb153A9adm8caMFZz+4Hv8Z8oyv1m2HHgycs65fVAnrSY3nh7dLNuheX1ueHkqlz/jN8vuL09GzjlXBgU3y955Tjc+/Ty6WfYZv1m2zDwZOedcGdWsIS5PuFn2rv/OZtBfP2K+3yy7z+Lq6fV+SXMlTZf0qqT0UN5c0lhJmyU9WmievpJmSMqV9LBCr3ySmkkaLWl++Ns0lCtMlxuW0yehriFh+vmh+3PnnCuzgptlH7zoCD5fvYVvPfw+D7w9jw3bdsUdWqUR15HRaKCHmfUEPgNuD+XbgV8CPy1inseBK4m6Cc8CBoTy24B3zSwLeDcMA5yVMO2VYX4kNQPuBI4GjgLuLEhgzjlXVpI4v3fbL2+WfWRMLifcN4Y/jprH2i074w4v5cWSjMzsbTMruHNsAtA2lG8xsw+IktKXJLUGGpvZxxY1W3kWOC+MHggMC++HFSp/1iITgPRQz5nAaDNba2briBJjQWJzzrn9UnCz7BvXncCJXVrwl3G5nPD7Mdz75hzyN3lT8OKkwjWjHwJvlTJNJpCXMJwXygBamdkKgPC3ZcI8S4uYp7jyb5B0paQcSTn5+flJrIpzzkW6t2nCY4P78vYNJ3FGt1b8/f2FnPD7Mdw1chZfbNheegXVTIUlI0nvSJpZxGtgwjR3ALuBF0qrroiy0pqsFDdP0nWZ2RNmlm1m2RkZGaUszjnnvimrVSMeurg37958Cuce0YbnJyzmpD+M5Y5XZ7B0rTcHL5BWURWb2WkljQ8NB84GTrXS7xjLI5zKC9oCy8P7lZJam9mKcBpuVcI87YqYJw84pVD5uFKW75xz+6VTiwbcP+gIrjs1i8fHL+CfOXm8PHEp5/fO5Jp+nenYokHcIcYqrtZ0A4BbgXPNrNSfBuH02yZJx4RWdJcCr4XRI4GCFnFDCpVfGlrVHQNsCPWMAs6Q1DQ0XDgjlDnnXIVr16w+vzv/cMbfcgrfP6YDI6ctp/8D47hh+BRyV1XfJuGK4zEWknKBOsCaUDTBzK4K4xYBjYHawHrgDDObLSkbeAaoR3SN6VozM0nNgRFAe2AJMMjM1oak9ShR44StwOVmlhOW8UPg52HZvzWzp0uLOTs723JycvZ73Z1zLtGqTdt58v3PeX7CYrbt2sNZPQ5iaL8surVpHHdo5ULSJDPLLnU6f6ZScjwZOecq0totO/nHB58z7KNFbNqxm9MOa8V1p3amZ9v0uEPbL56MypknI+fcgbBh2y6e+XAR//jwczZs28XJXTK4tn9nsjs2izu0MvFkVM48GTnnDqRN23fx/IQlPPn+QtZs2cmxBzfn2lM7c+zBzQkPoKkUPBmVM09Gzrk4bN25mxc/WcIT7y1k1aYdZHdoytD+nTm5S0alSEqejMqZJyPnXJy279rDP3OW8vi4BSzfsJ0j2jZhaP8sTjusZUonJU9G5cyTkXMuFezcvZd/T87jsXELWLJ2K4e1bszQfp05q8dB1KiReknJk1E582TknEslu/fs5bWpy/nLuFwW5m+hc8uGDO3XmbN7tiatZio86S3iyaiceTJyzqWiPXuNN2es4NExucxbuYmOzevzk36dOb93JrVSICl5Mipnnoycc6ls717j7dkreWTMfGYt30hmej2uPuUQBmW3pU5azdji8mRUzjwZOecqAzNj3Lx8Hh4znylL1nNQ47r8+OSDueSo9tStdeCTkiejcubJyDlXmZgZH+au4eEx8/n087W0aFiHH53Yie8f04EGdSrsGdnf4MmonHkycs5VVp8sXMOjY3N5f/5qmtavxRUndOLS4zrSuG6tCl+2J6Ny5snIOVfZTV6yjkfH5DJm7ioa1U3j8uM68sMTOpFev3aFLdOTUTnzZOScqypmLtvAo2Ny+d+sL2hQuyY/OLYj/3diJ1o0rFPuy/JkVM48GTnnqpp5X2zi0bG5vD59OXXSajD46A5cedLBtGpct9yW4cmonHkycs5VVQvyN/PY2AX8Z+oyatYQF2W346pTDiEzvd5+151sMoqrp9f7Jc2VNF3Sq5LSQ/npkiZJmhH+9k+Yp28oz5X0cOg8D0nNJI2WND/8bRrKFabLDcvpk1DXkDD9/ND9uXPOVVuHZDTkge8ewdibT+E7fTIZPnEJp9w/ltv+NZ3Fa7YckBjiuj13NNDDzHoCnwG3h/LVwDlmdjhRF+LPJczzOHAlkBVeA0L5bcC7ZpYFvBuGAc5KmPbKMD+SmgF3AkcDRwF3FiQw55yrzto3r8+9F/Rk3M/6cclR7fn3lGX0f2A8v3l9doUvO5ZkZGZvm9nuMDgBaBvKp5jZ8lA+C6grqY6k1kBjM/vYovOKzwLnhekGAsPC+2GFyp+1yAQgPdRzJjDazNaa2TqixFiQ2JxzrtrLTK/HPQN78MEt/bj8uI60a1a/wpd54O58Kt4PgZeLKP8OMMXMdkjKBPISxuUBmeF9KzNbAWBmKyS1DOWZwNIi5imu/BskXUl0VEX79u33ZZ2cc67Sa9m4Lr84u9sBWVaFJSNJ7wAHFTHqDjN7LUxzB7AbeKHQvN2B3wNnFBQVUU9pLS+KmyfpuszsCeAJiBowlLI855xzZVRhycjMTitpfGg4cDZwqiU06ZPUFngVuNTMFoTiPMKpvKAtUHA6b6Wk1uGoqDWwKmGedkXMkwecUqh8XPJr5pxzrrzF1ZpuAHArcK6ZbU0oTwfeAG43sw8LysNpuE2Sjgmt6C4FXgujRxI1diD8TSy/NLSqOwbYEOoZBZwhqWlouHBGKHPOOReTuFrTPQo0AkZLmirpr6F8KNAZ+GUon5pwDehq4EkgF1gAvBXK7wNOlzQfOD0MA7wJLAzT/x34CYCZrQV+DUwMr3tCmXPOuZj4Ta9J8ptenXNu36X0Ta/OOedcIk9GzjnnYufJyDnnXOz8mlGSJOUDi/ejihZEjzuq7KrKeoCvSyqqKusBvi4FOphZRmkTeTI6QCTlJHMRL9VVlfUAX5dUVFXWA3xd9pWfpnPOORc7T0bOOedi58nowHki7gDKSVVZD/B1SUVVZT3A12Wf+DUj55xzsfMjI+ecc7HzZOSccy52nowqmKQBkuZJypV0W+lzpCZJ/5C0StLMuGPZX5LaSRoraY6kWZKujzumspBUV9KnkqaF9bg77pj2l6SakqZIej3uWPaHpEWSZoSHPVfah1pKSpf0iqS54ftybIUty68ZVRxJNYHPiJ4mnkf0lPBLzKziO5QvZ5JOAjYTdeXeI+549kfo96q1mU2W1AiYBJxX2f4voTuVBma2WVIt4APgejObEHNoZSbpJiAbaGxmZ8cdT1lJWgRkm1mlvulV0jDgfTN7UlJtoL6Zra+IZfmRUcU6Csg1s4VmthMYDgyMOaYyMbP3gCrR1YaZrTCzyeH9JmAOxXQ9n8ossjkM1gqvSvvrMnSs+W2irmJczCQ1Bk4CngIws50VlYjAk1FFywSWJgznUQl3elWZpI5Ab+CTeCMpm3BaaypRD8ejzaxSrkfwEHALsDfuQMqBAW9LmiTpyriDKaODgXzg6XDq9ElJDSpqYZ6MKpaKKKu0v1yrGkkNgX8BN5jZxrjjKQsz22NmvYC2wFGSKuUpVElnA6vMbFLcsZST482sD3AWcE04zV3ZpAF9gMfNrDewBaiw696ejCpWHtAuYbgtsDymWFyCcI3lX8ALZvbvuOPZX+H0yThgQMyhlNXxwLnhWstwoL+k5+MNqezMbHn4uwp4leiUfWWTB+QlHG2/QpScKoQno4o1EciS1Clc/LsYGBlzTNVeuPD/FDDHzP4UdzxlJSlDUnp4Xw84DZgbb1RlY2a3m1lbM+tI9D0ZY2bfjzmsMpHUIDSMIZzWOgOodK1QzewLYKmkQ0PRqUCFNfJJq6iKHZjZbklDgVFATeAfZjYr5rDKRNJLwClAC0l5wJ1m9lS8UZXZ8cAPgBnhegvAz83szRhjKovWwLDQarMGMMLMKnWT6CqiFfBq9JuHNOBFM/tfvCGV2bXAC+HH9ELg8opakDftds45Fzs/Teeccy52noycc87FzpORc8652Hkycs45FztPRs4552Lnyci5FCbpIEnDJS2QNFvSm5K67GMdP6+o+JwrL96027kUFW7O/QgYZmZ/DWW9gEZm9v4+1LPZzBpWUJjOlQs/MnIudfUDdhUkIgAzmwp8IOl+STNDnzkXQdQ1hqT3Qh86MyWdKOk+oF4oeyGm9XCuVP4EBudSVw+ivpYKuwDoBRwBtAAmSnoP+B4wysx+G57KUN/M3pc0NDxM1bmU5cnIucrnBOAlM9sDrJQ0HjiS6FmI/wgPgf1POIpyrlLw03TOpa5ZQN8iyovqmqSgA8STgGXAc5IurcDYnCtXnoycS11jgDqSflRQIOlIYB1wUehYL4MoAX0qqQNRn0B/J3oqecHj/neFoyXnUpafpnMuRZmZSTofeEjSbcB2YBFwA9AQmEbUWeMtZvaFpCHAzyTtAjYDBUdGTwDTJU02s8EHej2cS4Y37XbOORc7P03nnHMudp6MnHPOxc6TkXPOudh5MnLOORc7T0bOOedi58nIOedc7DwZOeeci93/A4523kEycc02AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha = 0.001\n",
    "init_lambd = 0.1\n",
    "iterations = 1000\n",
    "lambd = [init_lambd/(i+1) for i in range(iterations)] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"stochastic\",\"train\")\n",
    "plt.title(\"stochastic gradient descent with decaying learning rate\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Based off my results with a decaying learning rate, it does seem like a better approach than a constant rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Kaggle submission"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Kaggle username: danish123, Kaggle best score: 0.91946\n",
    "#### Process to achieving my score: I just did hyperparameter tuning of learning rate, regularization amount, and iterations amount by training on the same training set for each combo and finding the combo with the best average validation accuracy. I saved the combos and their associated validation accuracies to a plain text file so I'd have my progress saved in case I encountered some error while the cell was running. I then did the straightforward process of training with all my training data with the best combo and producing test predictions with the optimized weights/theta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/cs189/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: divide by zero encountered in log\n",
      "/anaconda3/envs/cs189/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in add\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best validation accuracy:  0.7433333333333333\n",
      "best hyperparameters combo: (0.001, [0.01], 20000)\n",
      "CPU times: user 20min 32s, sys: 4.37 s, total: 20min 36s\n",
      "Wall time: 3min 27s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "alphas = np.geomspace(0.001,1.0,num=5).tolist()\n",
    "lambdas =  np.geomspace(0.01,1.0,num=5).tolist() #actually lambda\n",
    "iterations = [1000,10000,20000]\n",
    "theta = np.random.randn(n+1)\n",
    "best_hyperparams = None\n",
    "best_val_acc = 0.0\n",
    "\n",
    "for a in alphas:\n",
    "    for l in lambdas:\n",
    "        for i in iterations:\n",
    "            ll = [l]\n",
    "            theta, _ = \\\n",
    "            gradientdescent_reg_optimize(i,theta,ll,a,\"batch\",\"train\")\n",
    "            val_acc = \\\n",
    "            np.mean(predict(wine_xval,theta,\"val\") == wine_yval)\n",
    "            if val_acc > best_val_acc:\n",
    "                best_val_acc = val_acc\n",
    "                best_hyperparams = (a,ll,i)\n",
    "            file = open('hyperparams.txt','a')\n",
    "            file.write(','.join(map(str,(val_acc,a,l,i)))+\"\\n\")\n",
    "            file.close()\n",
    "\n",
    "print(\"best validation accuracy: \",best_val_acc)\n",
    "print(\"best hyperparameters combo:\",best_hyperparams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 413,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecXHXZ9/HPtT3Z9EJ6laAkNGGBoICIoAGVgKICKgmgERTbo7eAeCvWG+VRb7FhQKQIhiLNW3gwKEVvKQkhjZoQSjakbJJNNsn23ev54/wmObuZ2Z3ZndnZZL/v12tec+Z32jVnZs41vzJnzN0RERHJpoJ8ByAiIvsfJRcREck6JRcREck6JRcREck6JRcREck6JRcREck6JRfplJl908xuyHMMD5nZnB7e5yVmttHMdprZ8Bzva66Z/Sv22M3swDTXvcrM/pi76Pba33QzW5zlbV5nZv+ZzW3mipmdZGaVscfPm9lJ6SzbhX1l5biY2RlmtqC728lEUU/uTFIzs8eAP7p7Xk/iybj7jxLTZjYZeA0odvfmXOzPzK4CDnT3T8ViOC0X++oghmLgZ8BMd1/Wk/vOJTO7Cah09291YzPfB/5vbJuvA59x90e6ukF3v7gb8eSVu8/IxnbMbC7RcTw+tu2sHBd3f8DMfmRmh7n78mxsszOquUgbZpbTLxy53n4WjQLKgOczXdEi++Vny8zGAO8F7stgnR57zUMN8Kae2t8+5k/AvB7bm7vr1oUbMAG4B6gCtgC/CuUFwLeAN4BNwC3A4DCvDPhjWH4bsIjoJPZDoAWoB3YmtpVkn2cQney2AY8BB4fyy4G72y37C+DaMD0Y+D2wHlgH/AAoDPPmAv8L/BzYCvwgyX6vIqpVAbwJeIhzJ3BcKL8QeBGoBh4GJsXWd+ALwCrgtVh8a4Ea4FnghFA+C2gEmsL2l4Xyx4i+1XV2jCeH/c0JsW4GrozFcgywOOx3I/CzJM/3IGBX7Hn+I5S/K7xm28P9u2LrPBZex/8F6ohqXu23eznwKrADeAE4KzZvLvCvdsdsr22EeVOAx8N2FgK/Srw+Yf5dwIYQ5xPAjFA+LxzXxvC8/tJZXEn2fT7wSOzxrUBreM47gW/EXoOLwmvwREdxhXk3Ed57wElAJfC18PquBy5I83M5F7gpjeU6+8xcQPR+3gGsAT4XW+4kotpf4vHrwClhul94LtXhWP5Hu2WTHmvgYKLPf0s4jtvaH5fw+LPAaqLP6gPA2HbvmYuJPmfVwK8Bi81/N+Hz1xO3vJ+k98UbUAgsIzohlxMljePDvAvDiz8VGECUgG4N8z4H/AXoH7ZxFDAozHuMcPJMsc/ECe9UoDh8iFcDJcAkoDa2rcLwgZwZHt8H/C7EegDwTOLDEj6MzcAXiZpJ+yXZ91XsSS6Tw5u4KDb/zBDLwWEb3wL+3e5NvxAYltg+8ClgeFj+a0QnnbL2+4ttY/fx6eQYJ+K7nuiDfjjQwJ5E/CTw6TA9IHGMkjznNs8zxF4NfDrEfG54PDwW35vAjDC/OMk2PwaMJUqOnwiv55jY65BucnmSqMmuFDiR6EQVTy4XAgPD/P8Glsbm3US7LxAdxZVk39cAv25X9jrh5Nru2N1C9J7rl0lcRCfvZuB7RO/104ne30PT+GzOJb3k0tln5oPA2wAD3hOWPTIWX6rkcjXwz/B+mQCsbLds2u+BJMflZKIvS0eGY/hLQuKOvWf+BxgCTCT64jsrNn9YWGZQj5wne2In+9sNOC68cEVJ5v0d+Hzs8duJvi0WhQ/Xv4HDkqz3GB0nl/8E7ow9LiCqhZwUHv8LOD9Mnwq8GqZHEZ1c+8XWPRd4NEzPBd7s5PleRcfJ5SHgonax1RJqL2H5kzvZRzVwePv9JTs+nRzjRHzjY/OfAc4J008A3wVGdBJPm+dJlFSeabfMk8DcWHzfy/B9tBSYHXsdOk0u4aTRDJTHym5vf7xi84aEbSVqdjeRpHaaKq4k864Hrm5X9jrJk8vUDvaRMi6ik3ddu/fYJlJ8EWi33bmkkVw6+sykWPY+4Mux+FIllzW0PaHPiy+byXsgyXH5PfCT2LwB4X0/OfaeOT42/07g8tjj4rDMxEzep1297Zftwj1gAvCGJ+/QHkvUXJPwBtFJbxRRE8LDwAIze8vMfhI6jtPRZrvu3krUrDQuFN1OlDQAzguPIfqGVgysN7NtZraNqBZzQGzba9OMIZVJwC9i299K9I1vXGyZNvsws6+Z2Ytmtj2sMxgYkeb+OjrGCRti07VEH0SImmoOAl4ys0Vm9qEu7jOx35TPsT0zO9/MlsaO0yGk/5zjcVS7+652cST2UWhmV5vZq2ZWQ3Tio6P9ZBhXNVHtIx27j0cX4trS7vMVfw3bx/+bWOy/Ac5LPDazjjqvU31mMLPTzOwpM9satnt6B7HGjaXt+6DNe6ab74H254CdRE3s8fdgqvc97HndtqW5v25RcumatcDEFB2VbxGdbBMS3zQ3unuTu3/X3acTtd9/iKgNG6JvFB1ps10zM6Ikty4U3QWcZGbjgbPY80FZS1RzGeHuQ8JtkLcd4dLZvuOSLbuWqJltSOzWz93/nWw9MzsBuAz4OFFTxxCidnhLM56Ux7jT4N1Xufu5RMn1x8DdZlbe2XpJ9pnY77rY45Rxm9kkom/9lxI1pQ0hajKxVOuksB4Y2i7mibHp84DZwClECXtyIoRkMXYhruVEyTku1fOOl3cWV5e5++cT7zvg88DtsffhYR2smvQzY2alwJ+JRsSNCtt9MM1Y1xN9LhN2vzZpHOtMzwHlRE3L61Ku0dbBwOvuXpPm8t2i5NI1zxC9ia42s3IzKzOzd4d5fwK+amZTzGwA8CPgDndvNrP3mtmhZlZI1KHcRNSBB9GJcWoH+7wT+KCZvS/Udr5GlDT+DeDuVURNM38g6rR7MZSvB/4G/NTMBplZgZm9zcze08XnXkXUgRuP9TrgCjObAWBmg83sYx1sYyBRMqgCiszs28Cg2PyNwOQORlylPMadBW9mnzKzkaHml/gG19LROsGDwEFmdp6ZFZnZJ4DpRG3c6SgnOnlUhTguIPrWmhF3f4NoQMJ3zazEzI4HPhxbZCDR+2ILUd/ej9ptov37LNO4FgJHmllZB9tMprO4elyqzwxRP2Yp0TFpNrPTgPenudk7iT4LQ0PS+mJsXmfHeiMw3sxKUmz7duACMzsiJMAfAU+7++tpxvYeoibsHqHk0gXu3kL0gT6QqBO3kqhzDuBGouavJ4h+D1LPnjfYaOBuosTyItGIn8SP334BnG1m1WZ2bZJ9vkzUCf5Lok69DwMfdvfG2GK3E30zvL3d6ucTfWBeIGrWuBsY04WnjrvXEkZFhar9THe/l6gWsCA0eawEOvpdysNEb/JXiKr59bRtSrgr3G8xsyVJ1u/oGHdmFvC8me0kOubnuHt9Zyu5+xaimubXiE6Q3wA+5O6b09mpu78A/JSon2YjcCjRyLKuOA84lqj58TtEHecJtxAd03VEr/dT7db9PTA9vHb3ZRqXu28E/kFUC0n4L+BbYZtfT7FqZ3Hly16fGXffAXyJKFFUEx3vB9Lc3neJnudrRF/qbo1tt7Nj/Q+i0aAbzGyv95W7/52o7/XPRF9u3wack2ZcEDUB/i6D5bvFQkePiEhazGw6cDNwjOsEsk8wsw8TjZL8eI/tU+8NERHJNjWLiYhI1vW65GLRRfjWheF6S83s9Ni8K8xstZm9bGYfiJXPCmWrzezy/EQuIiIJva5ZzKKLFu509//brnw60SihY4jGez/CniGRrxD9CKqS6LIc54bOMxERyYN95SKCEI1OWeDuDcBrZraaKNEArHb3NQAWXVZ6NtGIlJRGjBjhkydPzmG4IiL7n2effXazu4/sbLnemlwuNbPzicbzf83dq4l+hRofvljJnl+mrm1XfmxnO5g8eTKLF2f1LylERPZ7Ztb+ShVJ5aXPxcweMbOVSW6zgd8Sjd8+gmgs908TqyXZlHdQnmy/88xssZktrqqqysIzERGRZPJSc3H3U9JZzsyuZ88voCtpe1mF8USXQ6CD8vb7nQ/MB6ioqOhdnU0iIvuR3jhaLP7L8bOIfu0N0S9kzzGzUjObAkwjugzLImBauBRICdEvVtP9Na2IiORAb+xz+YmZHUHUtPU60X+g4O7Pm9mdRB31zcAXwmVYMLNLiS4pUgjc6O4Z/3ugiIhkT68bitxTKioqXB36IiKZMbNn3b2is+V6XbOYiIjs+5RcREQk63pjn4uISJ/U0uo0tbTS2NJKU3MrTS3R46aWVprDvOYWp7k1mteSpGz3fYvTkqSsubWVC949hWHlqf42JjuUXESkz2htdRpbWmlobqWhuYWGpmi6MfG4ue3jppZWmpqjdZp235zG5naPdyeD2OPELcn6baabW3fPb+2BLnAzmH3EWCUXEdn/uDtNLU5dUwv1TS3UNbbsmd5d1rp7OrFMfRoJYc90Kw1NbR83trRm7TkUFxrFhQW7byWFRnFR28dFhQUUFxoDi4soSSxbFJWVxNYtLmr3uNAoCdsqKoj2U1RoFBVE84oKCyguMAoL9uwjPi+xTmGBJS3rCUouIpJSU0srtQ0t7GpsZldDM7saW6htaGZnQzO1jS3hvpldDS2759c1NlPX1EJdU2ubxFDXJkm00tKFr+lFBUZpUQGlxYWUFBZQWlwQPS4qpKQomh5QWhSmC8OyBZQUFu5ets28+OPiAkp3bzPaXkmqZFBomPXMSXpfpeQish9qbG5lR30TO+qbqQn3O+qbqKnb8zheXtsYkkO7RNLYnP43/f4lhfQvKaJfSQH9i4soKymkX3EBw8tL6De0kLKiwlAWbiWFlO2eLkg6v19xWKakkLKiAooKNQZpX6HkItJLNbe0sr2uieraJrbXNVK9q4nq2ka21Taxra4xKq9tSpos6ps6TwoDS4sYWFbEwLJiyksLKS8tYuTAUspLiygvKaJ/aSEDSoroX1pEeUk0v7w0SiADSovoX1IY3ZcW0a+4sMeaW2TfoOQi0kNqG5vZvKORqp0NbNnZwOadjWze2cDWXY1sq42SRfy+pr455bYKC4wh/YoZ3L+YQWXFDCwrYtzQfgwKySJxP7CsaPf8gWXFDOoX3Q8oLVIykJxSchHphvqmFqp2NLChpp6NNfVs3tHAll1R0qja0ciWXQ1s3tnA5h2N1DW1JN3GwNIihpQXM7R/CUP6lzB5RDlD+hUzpH8JQ/sXM7Q8Kh/SLyxTXszA0iK1+UuvpuQikkRrq1Nd27g7aWzYHiWQTTX1bKipZ8P2qLy6tmmvdQsMhpWXMmJACSMGlDJpYn+GDyhlxIA9ZSMGlDJiYAnDyksoLSrMwzMUyS0lF+mTGppbWL+tnnXb6lhXXUdluF+3rZZ12+rYsL2eppa2o5nMYHh5KaMHlzJ+aD+OmjSU0YPKGDW4jFGDyhg1qJSRA0oZ0r9ETU7S5ym5yH6pobmFtVvrWLu1NpY46lhXHSWPTTsaiF+z1QxGDSxj3NB+vHPCUMYcWsboQWW7k8foQWWMHFhKsUYriaRFyUX2WTX1Tby5pZY3ttTyxtZdvLmllte3RPfra+rbJI/iQmPM4H6MG9KPE6eNZNzQaHrc0H6MH9Kf0YPLKClS4hDJFiUX6dVaW523ttexetPO3bdVm3by2uZdbN3V2GbZEQNKmDisPzOnDmfi8P5MGt6ficP6M25Ifw4YWEqBmqpEeoySi/Qajc2tvLJxByvXbWf5uu08v247qzbtpLZxzyirEQNKeNvIAXxgxmgm704g5Uwc3p8BpXo7i/QWefk0mtnHgKuAg4Fj3H1xbN4VwEVAC/Ald384lM8CfkH0b5M3uPvVoXwKsAAYBiwBPu3ubb/SSq9U29jMs29U89SaLTy1ZisrKrfvvvbTwLIiDhk7mE8cPYFpBwzkwAMGcOABA3J+sT0RyY58fdVbCXwE+F280MymA+cAM4CxwCNmdlCY/WvgVKASWGRmD7j7C8CPgZ+7+wIzu44oMf22Z56GZGrt1lr+9sJG/vb8Bp59o5rmVqeowDhs/GDmvnsyh44bzKHjBjNpeH/9jkNkH5aX5OLuLwLJTh6zgQXu3gC8ZmargWPCvNXuviastwCYbWYvAicD54VlbiaqESm59CKbauq597l13PvcOl7asAOAd4weyGdPnMpxU4dz1KShlKtJS2S/0ts+0eOAp2KPK0MZwNp25ccCw4Ft7t6cZPm9mNk8YB7AxIkTsxSyJNPa6jz68iZuf/pNHnulipZW550Th/CtDx7MqdNHMWl4eb5DFJEcyllyMbNHgNFJZl3p7venWi1JmZP875i9g+WTcvf5wHyAioqKHvhbnr6nvqmFe5as44Z/rWFN1S4OGFjKvBOn8tEjx3PgAQPyHZ6I9JCcJRd3P6ULq1UCE2KPxwNvhelk5ZuBIWZWFGov8eWlBzU0t7DgmbX88h+r2byzgUPHDeYX5xzB6YeO0Q8PRfqg3tYs9gBwu5n9jKhDfxrwDFENZVoYGbaOqNP/PHd3M3sUOJtoxNgcIFWtSHKgpdW597l1/HzhK6zbVsfMqcP45bnvZObUYeqQF+nD8jUU+Szgl8BI4K9mttTdP+Duz5vZncALQDPwBXdvCetcCjxMNBT5Rnd/PmzuMmCBmf0AeA74fQ8/nT5r2dpt/Of9K1leuZ1Dxw3m6o8eyvEHjlBSERHMvW92PVRUVPjixYs7X1D2sq22kZ88/DJ/euZNRg4o5coPHswZh49VUhHpA8zsWXev6Gy53tYsJr3cwhc2csU9y6mubeLCd0/hK6dMY2BZcb7DEpFeRslF0rK9ronv/eUF/rykkoPHDOKWC49l+thB+Q5LRHopJRfp1FNrtvDVO5ayaUcDXzz5QL548jRdQVhEOqTkIim1tjq/eWw1P1v4CpOHl/PnS97FEROG5DssEdkHKLlIUlt2NvDVO5fxxCtVnHH4WH70kUN11WERSZvOFrKXJW9W8/k/LmFrbSM/POsQzjtmokaCiUhGlFykjXuWVHL5PSsYNaiUey55F4eMG5zvkERkH6TkIkD0S/trHn6Z6x5/lZlTh/GbTx6l/04RkS5TchF21DfxlQVL+ftLm/jksRO56owZuh6YiHSLkksft357HRf8YRGrNu3k+7Nn8OnjJuc7JBHZDyi59GGvbNzBnBufYUd9MzddcDQnTBuZ75BEZD+h5NJHPb1mC5+9ZTGlxYXc8bmZzBirjnsRyR4llz7owRXr+cqCpYwf1o+bLziGCcP65zskEdnPKLn0MTf972t8939e4MiJQ7nh/AqGakSYiOSAkksf4e78/JFVXPv3VZw6fRS/PPedlBUX5jssEdlPKbn0Ae7Ojx58kev/+RpnHzWeH3/0MAoL9It7EcmdvPyYwcw+ZmbPm1mrmVXEyiebWZ2ZLQ2362LzjjKzFWa22syutXA9EjMbZmYLzWxVuB+aj+fUW7W2Ot+6byXX//M1zj9uEj9RYhGRHpCvX8qtBD4CPJFk3qvufkS4XRwr/y0wD5gWbrNC+eXA3919GvD38FiA5pZWvn7XMm57+k0+956pfPeMGRQosYhID8hLcnH3F9395XSXN7MxwCB3f9Kj/2W+BTgzzJ4N3Bymb46V92mNza18acFz3PPcOv7PqQdx+ax36OKTItJjeuM1PqaY2XNm9riZnRDKxgGVsWUqQxnAKHdfDxDuD0i1YTObZ2aLzWxxVVVVLmLvFeqbWvjcrYt5cMUGvvXBg/nS+6YpsYhIj8pZh76ZPQKMTjLrSne/P8Vq64GJ7r7FzI4C7jOzGUCyM6NnGpO7zwfmA1RUVGS8/r5gV0Mzn71lMU+u2cKPzjqU846dmO+QRKQPyllycfdTurBOA9AQpp81s1eBg4hqKuNji44H3grTG81sjLuvD81nm7oX+b5re10TF960iOferOZnHz+cs945vvOVRERyoFc1i5nZSDMrDNNTiTru14Tmrh1mNjOMEjsfSNR+HgDmhOk5sfI+ZeuuRj55w1Msr9zGr887UolFRPIqX0ORzzKzSuA44K9m9nCYdSKw3MyWAXcDF7v71jDvEuAGYDXwKvBQKL8aONXMVgGnhsd9yqaaes6Z/ySrNu5k/qcrOO3QMfkOSUT6OIsGX/U9FRUVvnjx4nyH0W3rttXxyeufYtOOBm6YU8G73jYi3yGJyH7MzJ5194rOlutVzWL7gmsefolv3bci32EA8PrmXXz8uifZsquRWy86VolFRHoNXf4lQy+8VcOWXY35DoNXNu7gkzc8TUur86fPztR/3YtIr6Lksg9auW47n/790xQXFnDHvJlMGzUw3yGJiLSh5NIF+eymevaNaub+4RkGlRVz22eOZfKI8vwFIyKSgpJLhvL5S/d/v7qZz9y8mAMGlnLbZ2cybki/vMUiItIRJZd9xKMvb+LiW59l0vD+/PGiYzlgUFm+QxIRSUnJpQs88yvPdMvCFzbyhduWMG3UAG696FiG6d8jRaSXU3LJUE83ij20Yj1f/NNzzBg7iFsuPJbB/Yt7OAIRkcwpufRif1n2Fl+5YylHTBjCHy44mkFlSiwism9QcumCnhgtds+SSr5+1zIqJg/jxrlHM6BUL5WI7Dt0xspQTwwWu3PRWi67ZznHTR3ODXMq6F+il0lE9i06a/Uydyx6k8v+vIITpo3g+vMrKCsuzHdIIiIZ07XFuiBXzWL3PlfJ5fes4D0HjVRiEZF9mpJLxnLTLvbgivV87c5lzJwynN99+iglFhHZpym59AL/eGkjX/rTcxw5cSg3zFGNRUT2fUouXZDNVrElb1ZzyR+XMH3sIG684GjKNSpMRPYD+fonymvM7CUzW25m95rZkNi8K8xstZm9bGYfiJXPCmWrzezyWPkUM3vazFaZ2R1mltOfr2dztNgbW3bxmZsXM3pwGX+Yq9+xiMj+I181l4XAIe5+GPAKcAWAmU0HzgFmALOA35hZoZkVAr8GTgOmA+eGZQF+DPzc3acB1cBFPfpMumhbbSNz/7AId+cPc49m+IDSfIckIpI1eUku7v43d28OD58Cxofp2cACd29w99eA1cAx4bba3de4eyOwAJht0SWKTwbuDuvfDJzZU8+jq1pbna/esZTK6lquP7+CqSMH5DskEZGs6g19LhcCD4XpccDa2LzKUJaqfDiwLZaoEuVJmdk8M1tsZourqqq6HLB3cyzybx5bzaMvV/HtD02nYvKwbm1LRKQ3ylnvsZk9AoxOMutKd78/LHMl0AzcllgtyfJO8iToHSyflLvPB+YDVFRUdClDdLfL5bk3q/nZwlc44/CxfGrmpG5uTUSkd8pZcnH3Uzqab2ZzgA8B7/M9VYFKYEJssfHAW2E6WflmYIiZFYXaS3z5XqehuYX/uHs5owaV8cOzDsnrH4+JiORSvkaLzQIuA85w99rYrAeAc8ys1MymANOAZ4BFwLQwMqyEqNP/gZCUHgXODuvPAe7vqeeRqV8/+iqrN+3kRx85lIEaGSYi+7F8/ajiV0ApsDB8e3/K3S929+fN7E7gBaLmsi+4ewuAmV0KPAwUAje6+/NhW5cBC8zsB8BzwO9zGXhXKxsbttcz/4lX+eBhY3jv2w/IblAiIr1MXpKLux/YwbwfAj9MUv4g8GCS8jVEo8l6tf9+5BVaWp3LPvCOfIciIpJzvWG02D4n08Fia7fWctezlXxq5iQmDu+fm6BERHoRJZcMWRfGi93079cxYN6JU7MfkIhIL6TkkmM19U3csWgtHzpsDGMG98t3OCIiPULJpQs8g0tX/nX5enY2NDP33VNyGJGISO+i5JKhTEeL3bOkkreNLOfw8YNzE5CISC+k5JJDa7fWsuj1aj5y5Hj9YFJE+hQlly5Id7TYQyvXA3DG4WNzGI2ISO+j5JKhTCogj75UxTtGD2TCMA0/FpG+RcklR2rqm1j0+lZO0q/xRaQPSvsX+mb2LmByfB13vyUHMfV66bSK/e+qzTS3Oie/Q8lFRPqetJKLmd0KvA1YCrSEYgf6XHJJ90eUT67ZQv+SQt45cUjnC4uI7GfSrblUANO9u/+S1Ycser2aIycOpbhQLY8i0veke+ZbSfI//uqTOsuxNfVNvLShhorJQ3soIhGR3iXdmssI4AUzewZoSBS6+xk5iao3S6NVbMkb1bjD0foLYxHpo9JNLlflMoj9zZI3qiksMI6YoP4WEemb0kou7v64mY0Cjg5Fz7j7ptyF1bt11vG08q0aDhw5gPLSfP0Xm4hIfqXV52JmHyf6u+GPAR8Hnjazszteq8PtXWNmL5nZcjO718yGhPLJZlZnZkvD7brYOkeZ2QozW21m11q4noqZDTOzhWa2KtzntKMjnbFiL7xVw/Sxg3IZhohIr5Zuh/6VwNHuPsfdzyf658f/7MZ+FwKHuPthwCvAFbF5r7r7EeF2caz8t8A8YFq4zQrllwN/d/dpwN/D47zZsrOBDTX1TB+j5CIifVe6yaWgXTPYlgzW3Yu7/83dm8PDp4DxHS1vZmOAQe7+ZBgOfQtwZpg9G7g5TN8cK8+LF9fvAFDNRUT6tHQTxP8zs4fNbK6ZzQX+SpL/s++iC4GHYo+nmNlzZva4mZ0QysYBlbFlKkMZwCh3Xw8Q7lP+JN7M5pnZYjNbXFVV1fWIO+h0eWH9dgAOVs1FRPqwdDv0/8PMPgq8m6jbYb6739vROmb2CMl/G3Olu98flrkSaAZuC/PWAxPdfYuZHQXcZ2YzSN7VkfEPOt19PjAfoKKioks/CO3s0vkvrt/B6EFlDCsv6crmRUT2C2kPZ3L3PwN/zmD5Uzqab2ZzgA8B70v88t/dGwi/o3H3Z83sVeAgoppKvOlsPPBWmN5oZmPcfX1oPsvrKLZXq3YybdSAfIYgIpJ3HTaLmdm/wv0OM6uJ3XaYWU1Xd2pms4DLgDPcvTZWPtLMCsP0VKKO+zWhuWuHmc0Mo8TOB+4Pqz0AzAnTc2LlOZOqyuPurKnaxZQR5bkOQUSkV+uw5uLux4f7gVne76+AUmBhaGZ6KowMOxH4npk1E10g82J33xrWuQS4CehH1EeT6Ke5GrjTzC4C3iQaLp0zHTWKVe1sYGdDM1OVXESkj0vOJdYwAAAU1klEQVT7qsju/unOytLl7gemKE/Z9Obui4FDkpRvAd7XlTiybU3VLgCmjlSzmIj0bemOFpsRf2BmRcBR2Q9n35DqwpWvbU4kF9VcRKRv66zP5Qoz2wEcFu9vATbSA30bvVFHg8XWVO2ktKiAsYP79VxAIiK9UIfJxd3/K/S3XOPug8JtoLsPd/crOlq3L0p05hcUpPeHYiIi+6t0f+dyRbhm1zSgLFb+RK4C681SjRZ7c2utRoqJiJB+h/5ngC8T/b5kKTATeBI4OXeh9U6p6iTuzrptdZwwbWSPxiMi0hul26H/ZaLL7b/h7u8F3gl04/op+59ttU3UNrYwbqj6W0RE0k0u9e5eD2Bmpe7+EvD23IXVuyUbLLZuWx0A44YouYiIpHv5l8rwnyv3Ef3wsZo9l1/pU1JdW6yyOkou41VzERFJu0P/rDB5lZk9CgwG/l/OotoHVVZHV7FRchERSSO5mFkBsNzdD4HoL49zHlUv50nGi63bVkd5SSGD+xXnISIRkd6l0z4Xd28FlpnZxB6Ip9dLNVpsXXUd44b26/SS/CIifUG6fS5jgOfN7BlgV6LQ3c/ISVT7oHXb6tSZLyISpJtcvpvTKPYxyUaLbdhez+EThvR8MCIivVC6HfqPm9kkYJq7P2Jm/YHC3IbWSyVp9WpsbmXLrkZGDSzbe6aISB+U1u9czOyzwN3A70LROKJhyQJs3tkAwAGDSvMciYhI75Dujyi/ALwbqAFw91XAAd3ZsZl938yWm9lSM/ubmY0N5WZm15rZ6jD/yNg6c8xsVbjNiZUfZWYrwjrXWo571ds3i23aEZLLQCUXERFIP7k0uHtj4kH4P5dU129M1zXufpi7HwH8D/DtUH4a0QUypwHzgN+GfQ4DvgMcCxwDfCdcTJOwzLzYerO6GVtKlqRdbFNNPQAHqFlMRARIP7k8bmbfBPqZ2anAXcBfurNjd6+JPSxnT7KaDdzikaeAIWY2BvgAsNDdt7p7NbAQmBXmDXL3Jz36F69bgDO7E1umdtdc1CwmIgKkP1rscuAiYAXwOeBBd7++uzs3sx8C5wPbgfeG4nHA2thilaGso/LKJOU9ZlNNPWYwvLykJ3crItJrpVtz+aK7X+/uH3P3s939ejP7cmcrmdkjZrYyyW02gLtf6e4TgNuASxOrJdmUd6E8WTzzzGyxmS2uqsreRZ037WhgeHkpRYXpHk4Rkf1bumfDOUnK5na2kruf4u6HJLm1/4vk24GPhulKYEJs3niii2R2VD4+SXmyeOa7e4W7V4wc2bX/XUk2VGDTjgZ15ouIxHSYXMzsXDP7CzDFzB6I3R4FtnRnx2Y2LfbwDOClMP0AcH4YNTYT2O7u64GHgfeb2dDQkf9+4OEwb4eZzQyjxM4H2ievnNq0o179LSIiMZ31ufwbWA+MAH4aK98BLO/mvq82s7cDrcAbwMWh/EHgdGA1UAtcAODuW83s+8CisNz33H1rmL4EuAnoBzwUbjnj7cYib6ppYPqYQbncpYjIPqXD5OLubxCd+I/L9o7d/aMpyp3odzXJ5t0I3JikfDFwSFYDTKF9q5i7s2VXI8MHqOYiIpLQYXIxs3+5+/FmtoO2neRGlAf6/Nf1mvpmWlpdI8VERGI6q7kcH+4H9kw4+4Z4lq3eFf22dJiSi4jIbho7m6H2o8W21kbJZaiSi4jIbkou3bS75tJfyUVEJEHJpQvig8W2qllMRGQvSi4Zan/hymo1i4mI7EXJpZu27mqipLCA8pK++d9pIiLJKLl0gcfGi1XvamRoeTE5/gsZEZF9ipJLhpKNFhuqznwRkTaUXLpp665GdeaLiLSj5NIF8dFiUbOYkouISJySS4aSNYvpNy4iIm0puXRDa6uzva6Jof2L8x2KiEivouTSBYlWsR0NzbjDoH5KLiIicUouGdvTLlZT1wQouYiItKfk0g019SG5lCm5iIjE5SW5mNn3zWy5mS01s7+Z2dhQfpKZbQ/lS83s27F1ZpnZy2a22swuj5VPMbOnzWyVmd1hZjnvXU+MFqupawZgUL/O/tBTRKRvyVfN5Rp3P8zdjwD+B/h2bN4/3f2IcPsegJkVAr8GTgOmA+ea2fSw/I+Bn7v7NKAauCiXgcdHi6nmIiKSXF6Si7vXxB6W0/b/t5I5Bljt7mvcvRFYAMy26JorJwN3h+VuBs7MdrypJPpcBqvPRUSkjbz1uZjZD81sLfBJ2tZcjjOzZWb2kJnNCGXjgLWxZSpD2XBgm7s3tytPtc95ZrbYzBZXVVV1+znU1IdmMdVcRETayFlyMbNHzGxlkttsAHe/0t0nALcBl4bVlgCT3P1w4JfAfYnNJdmFd1CelLvPd/cKd68YOXJkV5/a7l0kai4DytTnIiISl7OzorufkuaitwN/Bb4Tby5z9wfN7DdmNoKoRjIhts544C1gMzDEzIpC7SVRnjPxbFZT38TA0iIKC3RFZBGRuHyNFpsWe3gG8FIoHx36UTCzY4ji2wIsAqaFkWElwDnAA+7uwKPA2WFbc4D7e+ZZRKPF9BsXEZG95as952ozezvQCrwBXBzKzwYuMbNmoA44JySQZjO7FHgYKARudPfnwzqXAQvM7AfAc8Dvcx387qHI9U0MVJOYiMhe8nJmdPePpij/FfCrFPMeBB5MUr6GaDRZj2gzFLmuSTUXEZEk9Av9bqipb9ZIMRGRJJRcuiAxHC2quahZTESkPSWXDFn8wpX1Taq5iIgkoeTSRa2tzs4GjRYTEUlGyaUL3J2djeG/XDRaTERkL0ouGUqMFqttaAGgvFTJRUSkPSWXLtrVGF1XrH9JYZ4jERHpfZRcusCJ1VxKVHMREWlPySVDibFiu2supaq5iIi0p+TSRbUhuajmIiKyNyWXLnCHXaFZTH0uIiJ7U3LJULhoM3WNIblotJiIyF6UXLpo1+5mMdVcRETaU3LpAnenNlFzUZ+LiMhelFy6aFdDM8WFRkmRDqGISHs6M3ZRbWML/YrVJCYikkzek4uZfd3M3MxGhMdmZtea2WozW25mR8aWnWNmq8JtTqz8KDNbEda5NvFXybm0q6FZl34REUkhr8nFzCYApwJvxopPA6aF2zzgt2HZYcB3gGOJ/nnyO2Y2NKzz27BsYr1ZuYzbiWouGoYsIpJcvmsuPwe+wZ7/3wKYDdzikaeAIWY2BvgAsNDdt7p7NbAQmBXmDXL3J93dgVuAM3MV8O4LVzaq5iIikkrekouZnQGsc/dl7WaNA9bGHleGso7KK5OUJ9vnPDNbbGaLq6qquhX/LtVcRERSyulXbzN7BBidZNaVwDeB9ydbLUmZd6F870L3+cB8gIqKiqTLpMWjmsuogWVd3oSIyP4sp8nF3U9JVm5mhwJTgGWh7308sMTMjiGqeUyILT4eeCuUn9Su/LFQPj7J8jmR+Jvj2oYW+g1XzUVEJJm8NIu5+wp3P8DdJ7v7ZKIEcaS7bwAeAM4Po8ZmAtvdfT3wMPB+MxsaOvLfDzwc5u0ws5lhlNj5wP25fg67Gpt10UoRkRR649nxQeB0YDVQC1wA4O5bzez7wKKw3PfcfWuYvgS4CegHPBRuOZP4Pxddbl9EJLlekVxC7SUx7cAXUix3I3BjkvLFwCG5ii/OLLr8S11zq2ouIiIp5Hso8j6psaWVllZXzUVEJAUlly5oaokGmqnmIiKSnJJLhuLjnvvpdy4iIkkpuXSDai4iIskpuXSD+lxERJJTcslQ/HrLqrmIiCSn5NINuraYiEhySi7doKsii4gkp+SSofj/kKnmIiKSnJJLNyi5iIgkp+TSDf3VoS8ikpSSS4YSrWJlxQUUFiT7KxkREVFyyVBhyC4ahiwikpqSS4YStRVd+kVEJDUllwwVhJqLOvNFRFJTcslQouZSVqzkIiKSSl6Ti5l93czczEaExyeZ2XYzWxpu344tO8vMXjaz1WZ2eax8ipk9bWarzOwOMyvJZcy7k0uRkouISCp5Sy5mNgE4FXiz3ax/uvsR4fa9sGwh8GvgNGA6cK6ZTQ/L/xj4ubtPA6qBi3IZd6JZrLRYlT4RkVTyeYb8OfANor+k78wxwGp3X+PujcACYLZFP5c/Gbg7LHczcGYugk0oDEdMzWIiIqnlJbmY2RnAOndflmT2cWa2zMweMrMZoWwcsDa2TGUoGw5sc/fmduWp9jvPzBab2eKqqqouxZ6ouSi5iIiklrMfa5jZI8DoJLOuBL4JvD/JvCXAJHffaWanA/cB02j7B5AJ3kF5Uu4+H5gPUFFRkU6NaS+JPpfSIjWLiYikkrPk4u6nJCs3s0OBKcCycBHI8cASMzvG3TfE1n/QzH4TOvsrgQmxzYwH3gI2A0PMrCjUXhLlOZPIZmXqcxERSanHz5DuvsLdD3D3ye4+mShxHOnuG8xsdOhHwcyOCfFtARYB08LIsBLgHOABd3fgUeDssPk5wP25jL+huRXQaDERkY70tmuYnA1cYmbNQB1wTkggzWZ2KfAwUAjc6O7Ph3UuAxaY2Q+A54Df5zLA+qaQXNTnIiKSUt6TS6i9JKZ/BfwqxXIPAg8mKV9DNJqsR9Q2RWMHdPkXEZHU1HGQobrGFkCXfxER6YiSS4aKCqJDNrR/Ti8EICKyT8t7s9i+5tKTD6SwAGYdkmyUtYiIgJJLxoaVl3DlB6d3vqCISB+mZjEREck6JRcREck6JRcREck6JRcREck6JRcREck6JRcREck6JRcREck6JRcREck6iy463PeYWRXwRhdXH0H0XzK9jeLKjOLKjOLKzP4a1yR3H9nZQn02uXSHmS1294p8x9Ge4sqM4sqM4spMX49LzWIiIpJ1Si4iIpJ1Si5dMz/fAaSguDKjuDKjuDLTp+NSn4uIiGSdai4iIpJ1Si4iIpJ1Si4ZMLNZZvayma02s8t7YH8TzOxRM3vRzJ43sy+H8qvMbJ2ZLQ2302PrXBHie9nMPpCr2M3sdTNbEfa/OJQNM7OFZrYq3A8N5WZm14Z9LzezI2PbmROWX2Vmc7oZ09tjx2SpmdWY2VfydbzM7EYz22RmK2NlWTtGZnZUeA1Wh3WtG3FdY2YvhX3fa2ZDQvlkM6uLHbvrOtt/qufYxbiy9tqZ2RQzezrEdYeZpfVf5SniuiMW0+tmtrQnj5elPjfk/f21m7vrlsYNKAReBaYCJcAyYHqO9zkGODJMDwReAaYDVwFfT7L89BBXKTAlxFuYi9iB14ER7cp+Alwepi8HfhymTwceAgyYCTwdyocBa8L90DA9NIuv1wZgUr6OF3AicCSwMhfHCHgGOC6s8xBwWjfiej9QFKZ/HItrcny5dttJuv9Uz7GLcWXttQPuBM4J09cBl3Q1rnbzfwp8uyePF6nPDXl/fyVuqrmk7xhgtbuvcfdGYAEwO5c7dPf17r4kTO8AXgTGdbDKbGCBuze4+2vA6hB3T8U+G7g5TN8MnBkrv8UjTwFDzGwM8AFgobtvdfdqYCEwK0uxvA941d07ugpDTo+Xuz8BbE2yz24fozBvkLs/6dGZ4JbYtjKOy93/5u7N4eFTwPiOttHJ/lM9x4zj6kBGr1341n0ycHc24wrb/Tjwp462ke3j1cG5Ie/vrwQll/SNA9bGHlfS8Yk+q8xsMvBO4OlQdGmo3t4Yq0anijEXsTvwNzN71szmhbJR7r4eojc/cEAe4ko4h7Yf+Hwfr4RsHaNxYToXMV5I9E01YYqZPWdmj5vZCbF4U+0/1XPsqmy8dsOBbbEEmq3jdQKw0d1Xxcp69Hi1Ozf0mveXkkv6krU39sg4bjMbAPwZ+Iq71wC/Bd4GHAGsJ6qWdxRjLmJ/t7sfCZwGfMHMTuxg2Z6Mi9CWfgZwVyjqDcerM5nGkqtjdyXQDNwWitYDE939ncD/AW43s0G52n8S2XrtchXvubT9EtOjxyvJuSHloin2n7PjpeSSvkpgQuzxeOCtXO/UzIqJ3jy3ufs9AO6+0d1b3L0VuJ6oKaCjGLMeu7u/Fe43AfeGGDaG6nSiGWBTT8cVnAYscfeNIca8H6+YbB2jSto2XXU7xtCZ+yHgk6EphNDstCVMP0vUn3FQJ/tP9RwzlsXXbjNRU1BRkni7JGzrI8AdsXh77HglOzd0sK2ef39l0kHTl29AEVFn1xT2dBTOyPE+jait87/blY+JTX+VqO0ZYAZtOznXEHVwZjV2oBwYGJv+N1FfyTW07Uz8SZj+IG07E58J5cOA14g6EoeG6WFZOG4LgAt6w/GiXQdvNo8RsCgsm+hwPb0bcc0CXgBGtltuJFAYpqcC6zrbf6rn2MW4svbaEdVk4x36n+9qXLFj9ng+jhepzw294v3l7kouGR2saMTFK0TfRq7sgf0dT1QVXQ4sDbfTgVuBFaH8gXYfwCtDfC8TG92RzdjDh2ZZuD2f2B5Ru/bfgVXhPvEmNeDXYd8rgIrYti4k6oxdTSwhdCO2/sAWYHCsLC/Hi6i5ZD3QRPRN8KJsHiOgAlgZ1vkV4YobXYxrNVHbe+J9dl1Y9qPhNV4GLAE+3Nn+Uz3HLsaVtdcuvG+fCc/1LqC0q3GF8puAi9st2yPHi9Tnhry/vxI3Xf5FRESyTn0uIiKSdUouIiKSdUouIiKSdUouIiKSdUouIiKSdUouIj3IzEab2QIze9XMXjCzB83soAy38c1cxSeSLRqKLNJDwkUO/w3c7O7XhbIjiH6Q+s8MtrPT3QfkKEyRrFDNRaTnvBdoSiQWAHdfCvzLov9TWRn+P+MTEF2+w8yeCP8LstLMTjCzq4F+oey2FPsRybuizhcRkSw5BHg2SflHiC7MeDgwAlhkZk8A5wEPu/sPzawQ6O/u/zSzS939iB6LWqQLlFxE8u944E/u3kJ04cHHgaOJru10Y7hA4X2hliOyT1CzmEjPeR44Kkl50r+P9ehPqk4kuvjhrWZ2fg5jE8kqJReRnvMPoNTMPpsoMLOjgWrgE2ZWaGYjiRLKM2Y2Cdjk7tcDvyf6q12AplCbEem11Cwm0kPc3c3sLOC/zexyoB54HfgKMIDoSroOfMPdN4T/V/kPM2sCdgKJmst8YLmZLXH3T/b08xBJh4Yii4hI1qlZTEREsk7JRUREsk7JRUREsk7JRUREsk7JRUREsk7JRUREsk7JRUREsu7/AxlSgV7GcKq4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = loadmat(\"data.mat\")\n",
    "wine_xfull, wine_yfull = data[\"X\"], data[\"y\"]\n",
    "m3, _ = wine_xfull.shape\n",
    "wine_xfull = np.concatenate([np.ones((m3, 1)), wine_xfull], axis=1)\n",
    "\n",
    "alpha, lambd, iterations = best_hyperparams\n",
    "\n",
    "theta = np.random.randn(n+1)\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"batch\",\"test\")\n",
    "plt.title(\"cost over iterations for all data (train + validation)\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 414,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test predictions saved\n"
     ]
    }
   ],
   "source": [
    "wine_xtest = data[\"X_test\"]\n",
    "m4, n4 = wine_xtest.shape\n",
    "wine_xtest = np.concatenate([np.ones((m4, 1)), wine_xtest], axis=1)\n",
    "                            \n",
    "predict(wine_xtest,theta,\"test\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# APPENDIX (don't run this code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.special import *\n",
    "from save_csv import *\n",
    "\n",
    "np.random.seed(1)\n",
    "\n",
    "def shuffle_train_val_split(val_amt=0, percent=0):\n",
    "    data = loadmat(\"data.mat\")\n",
    "    total_num_ex = data[\"X\"].shape[0]\n",
    "    shuffle = np.random.permutation(total_num_ex) #shuffle before split\n",
    "    data_xtrain, data_ytrain = data[\"X\"][shuffle], \\\n",
    "    data[\"y\"][shuffle] #split\n",
    "    if not val_amt: #split by fixed number or percentage\n",
    "        val_amt = int(percent * total_num_ex)\n",
    "    data_xval, data_yval = data_xtrain[:val_amt,:], data_ytrain[:val_amt,:]\n",
    "    data_xtrain, data_ytrain = data_xtrain[val_amt:,:], data_ytrain[val_amt:,:]\n",
    "    return data_xtrain, data_ytrain, data_xval, data_yval\n",
    "\n",
    "def sigmoid(values):\n",
    "    return expit(values)\n",
    "\n",
    "def predict(X,w,mode=\"train\"):\n",
    "    predictions = sigmoid(np.dot(X,w))\n",
    "    if mode != \"train\": #validation or test\n",
    "        predictions = (predictions >= 0.5).astype(int)\n",
    "        \n",
    "    if mode == \"test\":\n",
    "        results_to_csv(predictions)\n",
    "        print(\"test predictions saved\")\n",
    "    else:\n",
    "        return predictions\n",
    "\n",
    "def logres_cost_reg(X,y,theta,lambd):\n",
    "    m = X.shape[0]\n",
    "    regularization = (lambd/(2*m)) * np.sum(np.square(theta[1:]))\n",
    "    predictions = predict(X,theta)\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    J = ((-1.0/m) * (np.sum(np.dot(-y, np.log(predictions)) + \\\n",
    "                xlog1py(1-y, -predictions)))) + regularization\n",
    "    return J\n",
    "\n",
    "def gradientdescent_reg_optimize(iterations,theta,lambd,alpha,GD_type=\"batch\",mode=\"train\"):\n",
    "    x, y = None, None\n",
    "    if mode == \"train\":\n",
    "        x, y = wine_xtrain, wine_ytrain\n",
    "    elif mode == \"test\":\n",
    "        x, y = wine_xfull, wine_yfull\n",
    "        \n",
    "    costs = [0] * iterations\n",
    "    if len(lambd) == 1:\n",
    "        lambd = [lambd[0]] * iterations\n",
    "    next_theta = theta.copy()\n",
    "    \n",
    "    for i in range(iterations):\n",
    "        ###\n",
    "        if GD_type == \"batch\":\n",
    "            new_theta = batchGD_reg_iteration(x,y,theta,lambd[i],alpha)\n",
    "        elif GD_type == \"stochastic\":\n",
    "            new_theta = stochasticGD_reg_iteration(x,y,theta,lambd[i],alpha)\n",
    "        costs[i] = logres_cost_reg(x,y,theta,lambd[i])\n",
    "        theta = new_theta\n",
    "        \n",
    "    return theta, costs\n",
    "\n",
    "wine_xtrain, wine_ytrain, wine_xval, wine_yval = \\\n",
    "shuffle_train_val_split(percent=0.2)\n",
    "m1, n = wine_xtrain.shape\n",
    "m2, _ = wine_xval.shape\n",
    "wine_xtrain = np.concatenate([np.ones((m1, 1)), wine_xtrain], axis=1)\n",
    "wine_xval = np.concatenate([np.ones((m2, 1)), wine_xval], axis=1)\n",
    "\n",
    "alpha = 0.001\n",
    "lambd = [0.1] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "iterations = 1000\n",
    "\n",
    "def batchGD_reg_iteration(X,y,theta,lambd,alpha):\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    gradients = \\\n",
    "    (-1/len(y)*np.dot(X.T,(y-sigmoid(np.dot(X,theta))))).reshape(theta.shape)\n",
    "    gradients[1:] += (lambd/len(y) * theta[1:])\n",
    "    new_theta = theta - alpha*gradients\n",
    "    return new_theta\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"batch\")\n",
    "plt.title(\"batch gradient descent\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);\n",
    "\n",
    "alpha = 0.05\n",
    "lambd = [0.01] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "iterations = 1000\n",
    "\n",
    "def stochasticGD_reg_iteration(X,y,theta,lambd,alpha):\n",
    "    new_theta = theta.copy()\n",
    "    y = np.squeeze(y,axis=1)\n",
    "    np.random.shuffle(X)\n",
    "    for i in range(len(y)):\n",
    "        gradients = np.dot(y[i]-sigmoid(np.dot(X[i,:],theta)),X[i,:])\n",
    "        gradients[1:] += (lambd/len(y) * theta[1:])\n",
    "        new_theta = theta - alpha*gradients\n",
    "    return new_theta\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"stochastic\",\"train\")\n",
    "plt.title(\"stochastic gradient descent\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs)\n",
    "\n",
    "alpha = 0.001\n",
    "init_lambd = 0.1\n",
    "iterations = 1000\n",
    "lambd = [init_lambd/(i+1) for i in range(iterations)] #actually lambda\n",
    "theta = np.random.randn(n+1)\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"stochastic\",\"train\")\n",
    "plt.title(\"stochastic gradient descent with decaying learning rate\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);\n",
    "\n",
    "%%time\n",
    "\n",
    "alphas = np.geomspace(0.001,1.0,num=5).tolist()\n",
    "lambdas =  np.geomspace(0.01,1.0,num=5).tolist() #actually lambda\n",
    "iterations = [1000,10000,20000]\n",
    "theta = np.random.randn(n+1)\n",
    "best_hyperparams = None\n",
    "best_val_acc = 0.0\n",
    "\n",
    "for a in alphas:\n",
    "    for l in lambdas:\n",
    "        for i in iterations:\n",
    "            ll = [l]\n",
    "            theta, _ = \\\n",
    "            gradientdescent_reg_optimize(i,theta,ll,a,\"batch\",\"train\")\n",
    "            val_acc = \\\n",
    "            np.mean(predict(wine_xval,theta,\"val\") == wine_yval)\n",
    "            if val_acc > best_val_acc:\n",
    "                best_val_acc = val_acc\n",
    "                best_hyperparams = (a,ll,i)\n",
    "            file = open('hyperparams.txt','a')\n",
    "            file.write(','.join(map(str,(val_acc,a,l,i)))+\"\\n\")\n",
    "            file.close()\n",
    "\n",
    "print(\"best validation accuracy: \",best_val_acc)\n",
    "print(\"best hyperparameters combo:\",best_hyperparams)\n",
    "\n",
    "data = loadmat(\"data.mat\")\n",
    "wine_xfull, wine_yfull = data[\"X\"], data[\"y\"]\n",
    "m3, _ = wine_xfull.shape\n",
    "wine_xfull = np.concatenate([np.ones((m3, 1)), wine_xfull], axis=1)\n",
    "\n",
    "alpha, lambd, iterations = best_hyperparams\n",
    "\n",
    "theta = np.random.randn(n+1)\n",
    "\n",
    "theta, costs = \\\n",
    "gradientdescent_reg_optimize(iterations,theta,lambd,alpha,\"batch\",\"test\")\n",
    "plt.title(\"cost over iterations for all data (train + validation)\")\n",
    "plt.xlabel(\"Cost\")\n",
    "plt.ylabel(\"iteration\")\n",
    "plt.plot([i for i in range(iterations)], costs);\n",
    "\n",
    "wine_xtest = data[\"X_test\"]\n",
    "m4, n4 = wine_xtest.shape\n",
    "wine_xtest = np.concatenate([np.ones((m4, 1)), wine_xtest], axis=1)\n",
    "                            \n",
    "predict(wine_xtest,theta,\"test\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
